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ABSTRACT 


This  thesis  explores  the  possibility  and  feasibility  of  improving  existing  satel¬ 
lite  measurements  of  sea  surface  temperature  (SST)  by  the  incorporation  of  high- 
frequency  (HF)  radar-derived  surface  current  data.  Water  parcels  tagged  with  SST 
are  advected  using  particle  trajectories  calculated  by  integrating  surface  current  veloc¬ 
ity  data.  The  SST  of  these  advected  water  parcels  are  compared  to  SST  measurements 
at  the  final  times  and  locations  of  the  advected  water  parcels.  Different  methods  of 
generating  surface  currents  from  HF  radar  measurements  are  also  examined.  The 
Totals  current  method  is  a  local  fitting  method  which  generates  surface  current  mea¬ 
surements  by  solving  a  least-squares  equation  fitting  multiple  measurements  from 
different  radar  sites.  The  Open-boundary  Modal  Analysis  (OMA)  method  is  a  global 
method  which  fits  a  series  of  eigenfunction  modes  to  available  radial  measurements. 
These  modes  are  generated  by  solving  two  Laplacian  eigenvalue  problems  on  the  do¬ 
main  with  Dirichlet  and  Neumann  boundary  conditions,  and  adding  a  set  of  boundary 
modes  to  account  for  flow  across  open  boundaries.  Any  current  field  in  the  domain 
can  be  described  using  a  combination  of  these  modes.  The  two  methods  are  compared 
for  accuracy  against  an  analytic  solution  to  the  linear  Stommel  problem. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


vi 


TABLE  OF  CONTENTS 


I.  BACKGROUND .  1 

A.  REMOTE  SENSING .  1 

B.  CLOUDINESS  ON  THE  CENTRAL  CALIFORNIA  COAST  .  .  2 

C.  INCLUSION  OF  CURRENT  DATA  IN  SST  MODEL .  3 

D.  THESIS  OBJECTIVES  .  4 

II.  THEORY .  7 

A.  SST  MODEL .  7 

B.  DESCRIPTION  OF  HF  RADAR .  8 

C.  TOTALS  SURFACE  CURRENT  METHOD .  10 

D.  OPEN  MODAL  ANALYSIS  METHOD .  15 

1.  Calculation  of  Modes .  16 

2.  Unstructered  Grid .  21 

3.  Fitting  Current  Data  to  OMA  Modes .  21 

E.  PARTICLE  TRAJECTORIES .  24 

III.  METHODS .  27 

A.  MODEL  VALIDATION  AND  SST  COMPARISONS  .  27 

B.  DATA  USED .  30 

1.  SST  Data .  30 

2.  HF  Radar  Data .  31 

C.  CASE  STUDIES .  31 

1.  CENC  Domain  .  31 

2.  Dates  of  Study  .  31 

3.  Stommel  Ocean  Model .  32 

D.  MATLAB  PROCESSING  .  36 

1.  Processing  of  Radial  Data .  37 

2.  Processing  to  Totals  Currents .  37 

vii 


3.  OMA  Modes  Processing  .  38 

4.  Fitting  Radial  Data  to  OMA  Modes .  41 

5.  Particle  Trajectories .  41 

IV.  RESULTS .  45 

A.  STOMMEL  MODEL  CURRENTS  .  45 

1.  Current  Statistics .  45 

2.  Totals  Currents  Reconstruction  .  46 

3.  OMA  Currents  Reconstruction .  47 

4.  Mode  Coefficient  Penalty .  49 

B.  SST  COMPARISONS  AND  MODEL  TESTING .  51 

1.  Comparison  of  OMA  and  Totals  Advected  Methods  ...  57 

2.  Calculations  on  the  Distance  of  Advected  Particles  ....  60 

C.  RECONSTRUCTION  OF  SST  FIELD .  60 

D.  RECOMMENDATIONS  FOR  FURTHER  STUDY .  63 

1.  Refinement  of  Advection  Model .  63 

2.  Frontal  Advection .  63 

3.  Lagrangian  Coherent  Structures .  63 

4.  Adaptive  Interpolation  Functions  for  Scattered  Data  In¬ 
terpolation  .  64 

5.  Improvements  to  Current  Fits  from  OMA  Method  ....  66 

APPENDIX .  67 

LIST  OF  REFERENCES  .  69 

INITIAL  DISTRIBUTION  LIST  .  71 


viii 


LIST  OF  FIGURES 


1.  Cloudiness  on  the  central  California  coast .  2 

2.  SST  and  current  field  in  Monterey  Bay  on  19:10  GMT,  January  9,  2007.  4 

3.  Illustration  of  Totals  method .  9 

4.  GDOP  error  illustration .  11 

5.  Current  field  for  CENC  region  on  06:00  GMT,  19  October,  2006 .  13 

6.  GDOP  error  as  a  function  of  angle  separation .  15 

7.  Rectangular  and  Triangular  Grids .  22 

8.  Flow  chart  of  SST  comparison  procedure .  28 

9.  Advected  SST  comparison .  29 

10.  Static  SST  comparison .  29 

11.  Random  SST  comparison  .  30 

12.  Radar  coverage  for  the  CENC  region .  33 

13.  Stommel  Ocean  model  and  sample  HF  Radar  domain .  35 

14.  Triangular  Mesh  for  CENC  region .  39 

15.  Close-up  of  Triangular  Mesh  for  the  Monterey  Bay .  40 

16.  Dirichlet  modes  for  CENC  domain  .  42 

17.  Neumann  modes  for  CENC  domain .  43 

18.  Boundary  modes  for  CENC  domain .  44 

19.  Reconstruction  of  Stommel  current  held  using  Totals  method .  47 

20.  Reconstruction  of  Stommel  current  held  using  OMA  method .  48 

21.  Stommel  currents  reconstructed  from  different  OMA  modes  .  50 

22.  k  study  for  Stommel  domain .  51 

23.  SST  comparison  for  Totals  currents  for  January,  2007  .  52 

24.  SST  comparison  for  OMA  currents  for  January,  2007  .  52 

25.  SST  comparison  for  Totals  currents  for  October,  2006  .  53 

26.  SST  comparison  for  OMA  currents  for  October,  2006  .  53 


IX 


27.  SST  anomaly  comparison  for  Totals  currents  for  January,  2007  .  54 

28.  SST  anomaly  comparison  for  OMA  currents  for  January,  2007  .  54 

29.  SST  anomaly  comparison  for  Totals  currents  for  October,  2006  .  55 

30.  SST  anomaly  comparison  for  OMA  currents  for  October,  2006  .  55 

31.  SST  comparison  for  OMA  and  Totals  currents  for  January,  2007  ....  58 

32.  SST  anomaly  comparison  for  OMA  and  Totals  currents  for  January,  2007  58 

33.  SST  comparison  for  OMA  and  Totals  currents  for  October,  2006  ....  59 

34.  SST  anomaly  comparison  for  OMA  and  Totals  currents  for  October,  2006  59 

35.  Displacement  of  advected  particles  .  61 

36.  Recontructed  SST  field .  62 

37.  Lagrangian  Coherent  Structure .  64 

38.  Gaussian  weighting  function .  65 

39.  Gaussian  weighting  function,  slimmed  in  the  x  direction .  66 


x 


LIST  OF  TABLES 


I.  Statistics  for  Totals  Currents  calculated  on  Stommel  domain .  45 

II.  Statistics  for  OMA  Currents  calculated  on  Stommel  domain .  46 

III.  Satellite  Images  used  in  the  January,  2007  Case  Study  .  67 

IV.  Satellite  Images  used  in  the  October,  2006  Case  Study .  68 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


ACKNOWLEDGMENTS 


Thanks  to  my  advisors,  Frank  Giraldo  and  Jeff  Paduan  for  their  invaluable 
guidance.  Thanks  to  the  leadership  at  the  NOAA,  Southwest  Fisheries  Science  Center, 
Environmental  Research  Division,  for  making  my  studies  at  NPS  possible.  A  special 
thanks  also  to  Dr.  David  Kaplan  for  his  assistance  with  the  current  and  trajectory 
processing.  Without  David’s  MATLAB  toolboxes  and  help,  the  processing  in  this 
thesis  would  not  be  possible.  Most  of  all,  thanks  to  my  incredible  wife  Lisa,  who 
supported  me  during  the  long  nights  writing  this  thesis  by  feeding  me,  walking  the 
dog,  and  reminding  me  how  to  take  a  derivative.  Psalm  115:1. 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xiv 


I.  BACKGROUND 

A.  REMOTE  SENSING 

Remote  sensing  is  the  measurement  of  a  property  of  an  object  by  a  sensor  not 
in  direct  contact  with  the  object  being  measured.  Remote  sensing  techniques  extend 
important  tools  for  scientists  measuring  the  Earth’s  environment.  Many  established 
measurement  techniques  have  been  adapted  to  fit  on  remote  sensing  platforms,  in¬ 
cluding  satellites,  fixed  land-based  platforms,  Unmanned  Aerial  Vehicles  (UAVs),  Au¬ 
tonomous  Underwater  Vehicles  (AUVs),  or  the  underside  of  manned  aircraft.  These 
remote  sensing  techniques  offer  vast  improvements  over  traditional  measurement  and 
sampling  techniques  in  the  possible  spatial  domain  of  measurement  and  cost  effec¬ 
tiveness  of  measurement.  The  National  Oceanic  and  Atmospheric  Administration 
(NOAA)  has  launched  numerous  environmental  satellites  with  the  goal  of  providing 
remote  sensing  of  the  Earth  and  its  atmosphere. 

These  satellites  hold  multiple  sensors  onboard,  each  designed  to  measure  differ¬ 
ent  environmental  parameters.  Many  of  the  sensors  designed  to  take  a  measurement 
of  the  earth’s  surface  cannot  do  so  when  clouds  block  the  sensor’s  view  of  the  earth’s 
surface.  Passive  radiative  sensors,  such  as  the  Advanced  Very  High  Resolution  Ra¬ 
diometer  (AVHRR),  are  particularly  prone  to  this  problem.  The  AVHRR  sensor 
provides  a  measurement  of  sea  surface  temperature  (SST)  from  onboard  NOAA’s 
Polar  Operational  Environmental  Satellites  (POES).  Cloudiness  presents  significant 
problems  in  geographical  areas  with  persisting  cloud  cover  or  fog  banks,  such  as 
California’s  central  and  northern  coasts.  The  measurement  ability  of  satellites  in 
the  presence  of  clouds  might  be  improved  by  the  additional  information  provided  by 
measurements  that  are  not  affected  by  cloud  presence. 
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B.  CLOUDINESS  ON  THE  CENTRAL  CALIFORNIA  COAST 

To  illustrate  just  how  bad  the  cloud  problem  can  be  in  the  central  and  north¬ 
ern  California  coast,  statistics  were  calculated  on  the  region  covered  by  the  Califor¬ 
nia  Ocean  Currents  Monitoring  Program  (COCMP)  High  Frequency  Radars  (here¬ 
after  called  the  CENC  region).  Individual  satellite  pass  data  from  the  NOAA  POES 
AVHRR  sensor  were  downloaded  (courtesy  of  the  NOAA  Coast  Watch,  West  Coast 
Regional  Node)  and  masked  to  exclude  any  SST  pixels  not  within  the  CENC  region. 

The  percentage  of  cloudy  pixels  within  the  region  was  calculated  for  each  satellite 
pass. 


Cloudiness  Percentage 

(a) 


(b) 


Figure  1.  Cloudiness  on  the  central  California  coast  (CENC  region),  a.  Normal¬ 
ized  histogram  of  cloudiness  percentage  for  individual  satellite  passes,  b.  Average 
cloudiness  of  individual  satellite  passes  as  a  function  of  month. 
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A  total  of  1,966  satellite  images  were  analyzed  from  Jan  1,  2006  to  August 
10,  2007.  Figure  la  shows  a  histogram  (normalized  by  the  total  number  of  satellite 
images)  of  the  cloudiness  percentage  for  all  the  satellite  images  as  well  as  average 
cloudiness  as  a  function  of  month.  The  histogram  is  skewed  to  the  far  right,  indicating 
that  the  majority  of  SST  images  within  the  CENC  domain  have  vast  amounts  of  their 
data  destroyed  by  clouds.  Day  and  night  satellite  images  were  separated  in  Figure 
lb  due  to  the  different  cloudmasking  used  in  the  SST  processing.  The  nighttime 
cloudmasking  algorithm  only  uses  infrared  bands,  while  the  daytime  cloudmasking 
algorithm  utilizes  both  visible  and  infrared  bands.  Mean  cloudiness  per  image  is 
81.65%  for  the  year,  with  a  maximum  occurring  in  August  for  both  day  and  night 
images.  Upwelling  of  deeper,  colder  ocean  water  caused  by  seasonal  winds  often  reach 
a  peak  in  August.  The  colder  water  close  to  shore  provides  conditions  favorable  for 
fog  bank  formation. 

C.  INCLUSION  OF  CURRENT  DATA  IN  SST  MODEL 

As  part  of  the  California  Ocean  Currents  Monitoring  Program  (COCMP), 
multiple  land-based  high-frequency  (HF)  radar  sites  have  been  installed  along  the 
California  coast.  These  sites  form  a  sensor  network  designed  to  provide  a  continuous 
measurement  (in  space  and  time)  of  coastal  surface  ocean  currents.  The  presence  of 
clouds  do  not  affect  the  radar’s  ability  to  measure  ocean  currents.  We  hypothesize 
that  the  inclusion  of  surface  ocean  current  data  will  improve  existing  satellite  based 
SST  products  whose  measurements  are  hampered  by  the  presence  of  clouds  and  fog 
banks. 

The  surface  current  velocity  field  contains  information  relevant  to  the  SST 
field.  If,  for  example,  the  present  SST  field  is  unknown  due  to  clouds,  the  SST  field 
might  be  predicted  by  using  a  past  SST  field  obtained  from  some  satellite  pass  in  the 
past  when  the  sky  was  clear  and  the  surface  current  velocity  data  between  the  time 
of  the  past  measurement  and  the  current  time.  For  example,  if  there  was  a  strong 
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northward  current  running  along  the  coast,  we  might  expect  the  surface  temperature 
at  a  given  location  to  look  like  the  past  surface  temperature  at  some  point  to  the 
south.  Figure  2  shows  an  example  of  the  SST  and  current  fields  for  Monterey  Bay. 
Notice  the  stronger  westward  currents  off  Point  Pinos  (at  the  southern  boundary  of 
Monterey  Bay)  are  advecting  a  cold-water  filament  into  the  warmer  offshore  water 
mass. 
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Figure  2.  SST  (in  degrees  Celsius)  and  current  held  in  the  Monterey  Bay  on  19:10 
GMT,  January  9,  2007. 


D.  THESIS  OBJECTIVES 

The  main  objective  of  this  thesis  is  to  determine  if  the  inclusion  of  ocean 
surface  current  data  increases  the  ability  to  estimate  sea  surface  temperature  in  the 
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presence  of  clouds.  This  hypothesis  is  tested  by  building  a  model  of  SST  that  advects, 
or  moves,  parcels  of  surface  water  according  to  the  surface  currents  measured  by  the 
radars.  The  tracks  of  these  parcels  of  surface  water  are  called  particle  trajectories. 
Particle  trajectories  are  calculated  according  to  the  ordinary  differential  equation, 

dtts  .  /  \\  /T  \ 

—  =  u{t,x{t)),  (1.1) 

where  x(t)  is  the  position  vector  of  the  particle  composed  of  latitude  and  longitude 
components,  and  u(t,x(t))  is  the  velocity  vector  composed  of  latitudinal  and  longi¬ 
tudinal  velocities. 

In  the  particle  trajectories,  the  initial  time  and  location  of  a  particle  are  called 
the  departure  time  and  departure  point,  and  the  final  time  and  location  of  a  particle 
are  called  the  arrival  time  and  arrival  point.  In  the  model,  SST  at  the  arrival  point  at 
the  arrival  time  is  approximated  by  the  SST  of  a  water  parcel  which  has  been  advected 
from  the  departure  point  at  the  departure  time.  To  validate  the  model,  SST  at  the 
arrival  points  and  arrival  times  are  compared  against  SST  at  the  departure  points 
and  departure  times.  Note  that,  in  testing  the  model,  the  departure  and  arrival  times 
correspond  to  the  times  of  satellite  passes  to  allow  for  new  measurements  of  the  SST 
field.  Mathematically,  the  comparison  is  expressed  as: 

D  =  T(t0,  xq,  yo)  —  T(ti,xi,yi)  (1.2) 

where  T,  the  surface  temperature,  is  a  function  of  time  and  location  (in  latitude/longitude 
coordinates).  In  this  comparison  equation,  the  parcel  of  water  at  the  departure  point 
(x0,j/o)  at  the  departure  time  to  is  advected  by  the  surface  currents  to  the  arrival 
point  at  the  arrival  time  t\.  If  the  advected  water  parcels  temperature  is  a 

good  prediction  of  the  temperature  at  the  new  time  and  location,  then  this  difference 
equation  should  be  close  to  zero. 

As  a  secondary  thesis  objective,  the  different  methods  of  calculating  ocean 
currents  and  particle  trajectories  are  also  examined.  The  two  primary  methods  of 
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calculating  surface  currents  are  compared  to  each  other  and  to  an  analytic  solution 
to  the  Stommel  ocean  model  (discussed  in  chapter  III,  section  C). 

This  thesis  is  divided  into  four  chapters.  Chapter  I  presents  a  brief  background 
on  the  thesis  investigation.  Chapter  II  discusses  the  theory  behind  the  thesis  material. 
Chapter  III  discusses  the  methods  used  in  calculations  and  model  validation.  Chapter 
IV  discusses  the  results  of  the  model  validation  and  suggests  areas  of  further  research. 
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II. 


THEORY 


A.  SST  MODEL 

As  mentioned  above,  it  is  proposed  that  a  good  predictor  of  SST  is  a  model 
that  includes  both  past  data  about  the  SST  as  well  as  data  about  the  surface  current 
held.  In  reality,  the  surface  SST  held  is  much  more  complicated  and  will  depend 
on  additional  factors  such  as  surface  warming  or  cooling,  mixing  with  surrounding 
water  masses  (through  diffusion  or  turbulent  mixing),  and  horizontal  divergence  in 
response  to  three-dimensional  currents  (such  as  upwelling)  about  which  information  is 
not  known.  While  the  inhuence  of  these  processes  is  acknowledged,  they  are  ignored 
in  this  investigation.  The  model  being  investigated  is  the  one  in  which  horizontal 
advection  dominates  such  that  a  past  SST  held  and  surface  current  data  can  be  used 
to  estimate  a  present  SST  held. 

When  comparing  SST  measurements  to  validate  the  model,  both  absolute  SST 
and  SST  anomalies  are  compared.  SST  anomalies  are  dehned  as 


SSTa  =  T(tj,  Xi,  yi) 


( '52T(ti>xi’yi)\ 


2=1 

N 


(III) 


V  / 

where  N  is  the  total  number  of  SST  measurements  in  the  domain.  That  is,  the 
anomaly  is  the  difference  between  the  SST  measurement  at  some  point  and  the  average 
SST  of  the  region.  By  using  SST  anomalies,  it  is  hoped  that  some  of  the  inhuence 

of  large  scale  warming  or  cooling  is  removed.  That  is,  the  regional  change  in  time  of 

dT 

temperature  ( — due  tG  energy  hux  such  as  warming  by  the  sun  is  accounted 
for  by  only  looking  at  the  difference  between  the  local  point  temperature  and  the 
average  of  temperature  over  the  region.  If  the  entire  region’s  SST  is  changed  by  the 
same  amount  throughout,  then  the  SST  anomaly  at  any  point  within  the  domain  will 
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be  unchanged.  It  is  assumed  that  warming  or  cooling  will  occur  somewhat  uniformly 
over  the  region.  Localized  warming  or  cooling  will  produce  changes  in  the  anomalies. 


B.  DESCRIPTION  OF  HF  RADAR 

The  COCMP  seeks  to  monitor  ocean  currents  off  the  California  coast.  Ap¬ 
plications  of  surface  ocean  currents  include,  among  other  uses:  search  and  rescue 
operations,  prediction  of  chemical  and  oil  spill  trajectories,  and  prediction  of  trajec¬ 
tories  for  fish  larvae  transport. 

The  primary  instrument  currently  used  by  COCMP  is  the  Coastal  Ocean  Dy¬ 
namics  Application  (CODAR)  High  Frequency  Radar  antennae.  These  antennas  use 
high  frequency  electromagnetic  pulses  to  measure  ocean  currents  by  a  phenomenon 
known  as  Bragg  scattering  along  with  the  Doppler  effect  [1],  When  a  pulse  of  radar 
energy  is  sent  out,  it  scatters  off  the  ocean  surface  in  all  directions.  Some  energy  is 
backscattered  towards  the  radar  receiver  and,  of  that  energy,  some  is  Bragg  resonant. 
The  resonant  ocean  waves  are  those  traveling  in  a  direction  either  directly  towards  or 
directly  away  from  the  radar  antenna  with  a  wavelength  half  that  of  the  transmitted 
radiation.  For  these  ocean  waves,  the  antenna  picks  up  a  dramatic  spike  in  returned 
radiation.  This  is  due  to  two  reasons:  the  direction  of  travel  of  the  ocean  waves  opens 
up  more  of  the  ocean  wave  face  to  reflect  radiation  back  towards  the  radar,  and,  more 
importantly,  the  reflected  radiation  off  one  resonant  ocean  wave  is  in  phase  with  the 
reflected  radiation  off  the  next  resonant  ocean  wave.  This  phenomenon  is  known  as 
Bragg  scattering. 

From  the  wavelength  of  the  resonant  ocean  waves,  the  theoretical  speed  of  the 
resonant  ocean  wave  can  be  calculated  from  the  deep-water  dispersion  relationship. 
The  deep-water  dispersion  relationship  is: 


where  C  is  the  wave  speed,  g  is  the  acceleration  due  to  gravity,  and  A  is  the  wavelength 


of  the  ocean  wave  [2].  The  wavelength,  theoretical  speed,  and  the  direction  of  travel 
for  the  Bragg-resonant  ocean  waves  are  known.  The  received  signals  exhibit  a  Doppler 
shift  that  is  slightly  different  than  the  shift  expected  due  to  the  theoretical  speed  of  the 
ocean  waves.  In  the  absence  of  any  other  influences,  the  Doppler  shift  would  always  be 
equal  to  the  speed  of  the  traveling  wave.  The  difference  between  the  theoretical  wave 
speed  and  the  speed  measured  by  Doppler  shift  is  due  to  ocean  currents.  Estimates 
of  the  uncertainty  of  radial  speed  measurements  are  in  the  range  of  6. 9-9. 2  cm/sec 


[3], 


Because  the  radars  measure  the  ocean  current  speed  using  a  Doppler  shift,  a 
single  Radar  site  can  only  give  velocity  data  for  currents  moving  towards  or  away 
from  that  site.  A  single  radar  can  only  measure  the  component  of  current  which 
points  towards  or  away  from  the  radar.  Therefore,  measurements  from  multiple  sites 
are  needed  at  any  one  location  to  give  a  complete  velocity  picture.  Figure  3  shows 
an  illustration  of  how  a  Total  current  vector  is  made  from  the  individual  radial  com¬ 
ponents. 


Figure  3.  Illustration  of  Totals  method.  Figure  taken  from  [1], 
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The  angles  between  sites  should  be  as  close  to  orthogonal  as  possible  to  pre¬ 
vent  a  problem  known  as  the  geometric  dilution  of  precision  (GDOP).  GDOP  occurs 
when  the  angles  between  radar  sites  are  small.  Small  errors  in  the  radial  velocity 
measurements  can  lead  to  large  errors  in  the  final  current  measurement  [4],  Figure 
4  illustrates  the  problem  that  can  arise  when  sites  with  small  separation  angles  are 
used  to  make  a  total  current  vector  with  sampling  error.  When  the  measurements 
are  collected  without  any  error,  then  the  resulting  current  is  also  without  error.  But 
when  sampling  error  is  present  and  the  angle  spread  between  sampling  sites  is  small, 
even  small  errors  (10%  in  the  figure)  can  lead  to  large  errors  in  the  final  current. 
GDOP  is  discussed  further  in  later  sections. 

C.  TOTALS  SURFACE  CURRENT  METHOD 

The  processing  of  radar  measurements  to  surface  current  velocities  and  particle 
trajectories  follows  two  methods.  The  first  is  what  will  be  referred  to  as  the  ‘Totals 
method’,  and  the  second,  the  ‘Open-boundary  Modal  Analysis  (OMA)  method’.  Both 
methods  were  used  to  generate  currents  in  this  thesis.  The  resulting  currents  are 
compared  to  each  other  and  to  the  analytic  solution  of  the  Stommel  ocean  model. 

In  the  processing  of  radial  current  measurements  to  currents  by  the  Totals 
method,  a  regular  latitude/longitude  grid  is  specified  where  current  measurements 
are  desired.  All  radial  measurements  that  fall  within  a  specified  radius  of  the  grid 
location  are  grouped  together.  A  least  squares  fit  is  obtained  for  the  total  currents 
at  each  grid  location.  The  least  squares  problem  is  formed  as: 

r  =  Au.  (II. 3) 

If  a  grid  point  has  n  associated  radial  observations,  then  r  is  an  nxl  vector 
of  the  measured  radial  velocities,  u  is  a  column  vector  (size  2x1)  of  the  complete 
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Orthogonal  measurements,  no  error 


Orthogonal  measurements,  with  10%  error 


Close  measurements,  no  error 


/  \ 

I  \ 

i  \ 


Close  measurements,  with  10%  error 


/  \ 

/  \ 

/  \ 


Figure  4.  Example  of  GDOP  in  radial  measurements  of  a  current  vector  from  or¬ 
thogonal  measurement  sites  (top  panes)  and  measurement  sites  that  are  separated  by 
30°  (bottom  panes).  Black  vectors  are  the  radial  measurements  and  red  vectors  are 
the  resulting  total  fitted  current.  In  the  left  panes,  radial  measurements  are  sampled 
without  error,  and  in  the  right  panes  radial  measurements  are  sampled  with  a  10% 
error. 
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velocity  composed  of  u  and  v  components,  and  A  is  the  transformation  matrix 


-  cosdi 
cos  02 


sindi  ^ 
sin  02 


A 


(II.4) 


y  cos  6n  sindn  J 

where  9i  is  the  bearing  from  the  ith  radial  measurement  to  the  radar  site  that  took 
the  measurement  [5] .  At  a  minimum,  a  radial  measurement  from  two  different  radar 
sites  are  needed  to  generate  the  total  current  measurement  at  a  single  point.  Ra¬ 
dial  measurements  from  additional  radar  sites  will  provide  more  information  about 
the  currents  at  the  grid  point  and,  in  the  absence  of  problems,  give  a  better  mea¬ 
surement  of  the  total  current.  The  surface  current  field,  which  now  contains  current 
measurements  at  the  specified  points  of  the  regular  grid,  may  be  cleaned  by  remov¬ 
ing  currents  with  unreasonably  large  surface  current  magnitudes  or  with  large  errors. 
Figure  5  shows  an  example  one- hour  current  field  of  the  radar  domain  off  central 
California.  This  current  field  was  calculated  using  the  Totals  method. 

Error  measurements  are  also  generated  in  calculating  totals  currents.  Most 
error  measurements  center  around  the  concept  of  GDOP  and  are  a  function  of  the 
spread  of  the  angles  of  radial  measurements  that  make  the  total  current  vector.  If  A 
is  the  matrix  which  is  used  in  the  least  squares  fit  of  radial  measurements  to  total 
currents,  then  the  covariance  matrix  of  the  total  current  components  u  and  v  is  the 
2x2  matrix: 


Cov  =  cr2(ATA)  1, 


(II.5) 


where  a2  is  the  variance  of  the  radial  measurements.  The  variance  of  the  u  and  v 
total  velocity  components  are  the  first  and  second  diagonal  entries  of  this  matrix, 
respectively.  The  covariance  of  u  and  v  are  the  off-diagonal  entries  [6].  Notice  that 
as  the  angles  between  two  radar  sites  become  close  to  parallel,  the  matrix  (ATA) 
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Longitude 


Figure  5.  Totals  current  field  for  CENC  region  on  06:00  GMT,  19  October,  2006. 
Every  4th  vector  is  plotted. 
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becomes  close  to  singular,  and  the  inverse  of  ( AT A )  contains  very  large  values.  For 
angles  which  are  exactly  parallel,  the  inverse  of  (ATA)  does  not  exist.  As  the  angles 
between  two  radar  sites  become  close  to  parallel,  any  sampling  error  will  produce 
very  large  variances  in  the  u  and  v  fitted  velocities.  That  is,  even  if  a2  is  very  small, 
the  diagonal  elements  of  equation  II. 5  will  still  tend  towards  infinity  for  close  radar 
measurements. 

A  GDOP  error  estimate  is  formed  by  setting  a  unit  radial  measurement  vari¬ 
ance  (cr2  =  1).  With  this  substitution  for  a,  the  error  estimate  is  purely  a  measure 
of  the  loss  of  precision  due  to  the  angular  spread  of  the  radial  measurements.  If 
additional  information  is  known  about  the  variance  of  the  radial  measurements  that 
makes  the  total  current  (cr2),  then  that  can  be  added  to  equation  II. 5  to  produce 
another  type  of  error  measurement. 

The  covariance  matrix  is  formed  as  above  in  equation  II. 5,  and  the  total  GDOP 
error  measurement  is  defined  to  be: 

E  =  ^||  Cov  || .  (II. 6) 

That  is,  the  total  GDOP  error  is  the  square  root  of  the  covariance  matrix  2- 
norm,  or  the  square  root  of  the  largest  eigenvalue  of  the  covariance  matrix.  This  total 
GDOP  error  measurement  is  a  function  of  angle  separation.  Figure  6  is  a  plot  of  the 
GDOP  against  angle  separation.  The  GDOP  error  function  reaches  a  minimum  of 
one  at  orthogonal  angles  (angle  separations  of  90°  or  270°)  and  tends  towards  infinity 
as  the  angle  vectors  approach  parallel  (0°  or  180°).  This  total  error  is  close,  but  not 
equal,  to  the  square  root  of  the  sum  of  the  squares  of  the  two  variances  (diagonal 
elements). 

Another  error  measurement  involves  finding  the  GDOP  using  only  the  two 
most-orthogonal  radial  measurements  that  make  a  current  vector.  This  is  perhaps 
more  useful  than  the  GDOP  calculated  with  all  the  radial  measurement  angles.  In 
the  full  GDOP  error  measurement,  additional  radial  measurements  with  the  same 
angles  (additional  rows  in  the  A  matrix)  will  lower  the  GDOP  error.  The  reduction 
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Figure  6.  GDOP  error  as  a  function  of  angle  separation 

in  GDOP  is  somewhat  artificial  because  the  additional  radial  measurements  measure 
the  same  component  of  current  and  add  no  new  information  to  the  least-squares  fit. 

D.  OPEN  MODAL  ANALYSIS  METHOD 

The  second  method  of  processing  radar  measurements  to  the  current  veloc¬ 
ity  field  uses  a  method  called  Open-boundary  Modal  Analysis,  or  OMA.  The  OMA 
method  used  in  this  thesis  follows  the  procedure  described  in  [7].  The  general  idea 
of  OMA  is  to  generate  a  set  of  modes  for  a  given  domain  which  can  be  used  to 
approximate  any  current  field  on  that  domain.  The  modal  series  approximation  is 
determined  by  minimizing  a  cost  function  to  find  the  ideal  combination  of  the  modes 
which  gives  the  best  fit  to  available  measurement  data.  These  modes  depend  only  on 
the  shape  of  the  domain.  Once  they  are  calculated,  they  can  be  stored  for  repeated 
use  on  the  same  domain. 
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In  OMA,  the  domain  of  interest  need  not  be  rectangular.  Furthermore,  the 
boundary  of  the  domain  can  be  composed  of  multiple  open  and  closed  segments. 
An  open  segment  is  one  along  which  the  normal  component  of  the  velocity  can  be 
nonzero.  Vector  currents  are  allowed  to  flow  into  or  out  of  an  open  segment.  A  closed 
segment  is  one  in  which  the  normal  component  of  the  velocity  is  fixed  at  zero  and  no 
current  is  allowed  to  flow  through  the  segment,  although  currents  can  flow  parallel  to 
a  closed  segment. 

In  theory,  the  OMA  method  offers  several  advantages  over  the  totals  currents 
method.  The  generation  of  the  OMA  modes  is  only  done  once.  Afterwards,  when 
fitting  the  currents  to  the  OMA  modes,  one  least-squares  matrix  equation  needs  to 
be  solved.  This  is  computationally  similar  to  the  problem  of  finding  currents  with 
the  totals  method.  Unlike  the  totals  method  however,  the  modes  are  defined  over 
the  entire  domain,  even  for  sparse  current  measurements.  The  fit  results  in  a  current 
field  which  is  also  defined  over  the  entire  domain  without  gaps.  This  is  a  nice  feature 
of  the  OMA  method,  but  should  be  utilized  with  caution.  Currents  will  be  reported 
even  in  areas  with  few  or  no  actual  measurements.  Currents  in  these  areas  are  an 
extrapolation  of  the  modal  fit  to  data  in  other  areas,  and  they  do  not  represent  the 
“real”  currents.  Attention  should  be  paid  to  how  much  real  data  goes  into  making 
the  fit. 

1.  Calculation  of  Modes 

The  first  step  in  OMA  is  to  calculate  the  basis  functions  (modes)  that  are  later 
used  to  reconstruct  the  current  velocities.  The  modes  depend  only  on  the  shape  of  the 
domain,  and  the  calculation  of  the  modes  needs  to  be  done  only  once  for  a  particular 
domain.  Proceeding  as  in  [7],  it  is  noted  that  any  vector  field  u  on  a  domain  0  can 
be  decomposed  into  an  irrotational  (curl-free)  vector  component  uv  and  a  solenoidal 
(divergence-free)  vector  component  u^.  This  decomposition  is  known  as  Helmholtz’s 
Theorem,  or  the  Fundamental  Theorem  of  Vector  Analysis,  and  is  written  as: 

u  =  ulfi  +  u^.  (II. 7) 
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Because  is  curl-free,  or  irrotational,  it  is  a  conservative  field  and  there  exists  a 
scalar  potential  function  ip  such  that: 

u v  =  V</?.  (II. 8) 

The  field  is  divergence-free,  or  solenoidal,  and  can  be  expressed  in  terms  of  a 
stream  function  xp: 

=  k  x  Vip,  (II. 9) 

where  k  is  the  unit  vector  pointing  out  of  the  plane. 

It  can  be  verified  that  u v  is  curl-free  and  is  divergence-free.  By  taking  the 
divergence  of  u,  the  following  nonhomogeneous  partial  differential  equation  results 
for  Lp\ 

_  „  .  .  „  „  „  „  .  d  -  d  ,dp>~  dtp 

V  u  =  V  (uv  +  itip)  =  V-uv+V  u*  =  =  v  <P- 

(II-10) 

By  taking  the  curl  of  u,  the  following  nonhomogeneous  PDE  results: 

V  x  u  =  V  x  (uv  +  u^p)  =  Vx(h  Vt/0  =  V  x  +  ~^-j)  =  'V2'ipk-  (II. 11) 

Taking  the  dot  product  of  equation  II.  11  and  k  gives 

k  •  (V  x  u)  =  V2tp.  (11.12) 

a.  Interior  Modes 

This  formulation  for  p>  and  xp  is  not  unique.  For  example,  <p  and  ip 
plus  any  constant  still  satisfy  equations  II.  10  and  11.12.  Hence  there  is  some  degree 
of  freedom  and  some  flexibility  on  the  choice  of  boundary  conditions  along  dfl,  the 
boundary  of  the  domain  0.  The  authors  of  [7]  chose  to  define  xp  to  be  zero  along  the 
entire  boundary,  and  specify  the  normal  component  of  the  flow  of  u  solely  in  terms 
of  ip.  To  solve  this  pair  of  equations  (II.  10  and  11.12),  the  method  of  eigenfunction 
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expansion  is  used.  p  and  xp  are  expanded  in  terms  of  their  eigenfunctions,  hereafter 
called  modes,  which  are  solutions  to: 

V2^  =  A  iXpi,  xpildn  =  0  (11.13) 

V2(pi  =  \itpi,  (n  •  V<ft)|dQ  =  0.  (11.14) 

Because  of  the  specific  boundary  conditions  imposed  on  ip  and  xp,  the 
eigenmodes  of  p  are  called  Neumann  modes  and  the  eigenmodes  of  xp  are  called 
Dirichlet  modes.  Equations  11.13  and  11.14  are  two  Sturm- Liouville  equations,  and 
<pi  and  p>i  are  eigenfunctions  of  the  Sturm-Liouville  equations.  They  are  linearly 
independent  and  form  a  basis  for  the  space  of  solutions  of  the  Sturm-Liouville  problem. 
As  a  result,  the  vector  fields  'Vpi  and  k  x  VA  reconstructed  from  the  modes  p,.  and 
xpi  form  a  basis  for  all  two  dimensional  vector  fields  in  fl  [8] .  The  series  of  eigenmodes 
of  (p  and  xp  are  used  to  reconstruct  a  velocity  held  u  as: 

OO  OO 

u  =  J>f  (fc  x  VxPi)  +  5>fV^,  (11.15) 

1=1  1=1 

where  the  cq  are  unknown  coefficients. 

The  important  limitation  so  far  being  that  any  combination  of  the 
currents  reconstructed  from  <p  and  xp  will  have  zero  how  on  the  domain  boundary  dll. 
This  is  verihed  by  checking  the  normal  component  of  the  current  associated  with  each 
mode,  h  ■  u,  on  the  boundary. 


h  ■  =  n  ■  'Vpi  =  0, 


(11.16) 


due  to  the  boundary  conditions  of  equation  11.14.  Also, 

~  tt  t-,  ,  \  ~  9xp- 

n  u^  =  n-  (k  x  Vxpi)  =  n  ■ 


(11.17) 


To  verify  that  this  is  zero  on  the  domain  boundary,  let  t  be  a  unit  vector  tangent 
to  the  domain  boundary.  Because  xppdQ  =  0,  the  boundary  is  a  contour  line  in  xp , 
and  'Vxp  is  necessarily  orthogonal  to  the  boundary  and  t.  If  t  is  composed  of  vector 
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components,  (xti  +  ytj ) ,  then 

-  ,  .  -  .dip-  dtp-.  dip  dip  ,  . 

t  ■  Vipi  =  (xti  +  ytj)  ■  Vipi  =  (xti  +  ytj)  •  (— i  +  —j)  =  xt—  +  =  0,  (11.18) 


which  implies 


dip 

hw-  = 
ox 


-yt 


dip 
<  >!/  ’ 


(11.19) 


Combining  this  with  equation  11.17  and  the  specification  that  t  =  n  x  k,  gives 
h  =  (— yti  +  xtj )  (due  to  the  orthogonality  of  n  and  t).  Thus, 


/  ^  ,  dip*  dip^  dip  dip  dip  dip  .  . 

n  u *  =  (-y‘l+x‘3)^-ayl+di3)  =  ytay+Xt^  =  -x'Tx+Xtdx  =  °-  (IL20) 

And  so  it  is  shown  that  the  normal  component  of  the  current  associated  with  ipi  along 

the  domain  boundary  is  also  zero. 


b.  Boundary  Modes 

Since  an  arbitrary  domain  will  consist  of  both  closed  and  open  boundary 
segments,  the  full  reconstruction  of  a  current  field  must  include  flow  across  open 
boundary  segments.  Open  boundary  flow  is  accounted  for  by  including  a  set  of 
boundary  modes,  ip\.  These  boundary  modes  satisfy: 


VV |n  =  ^  <f  g(s)ds, 
A  JdQ 

(n  ■  \dQ  =  n  ■  ub\dn  =  g(s), 

JJydA= °- 


(11.21) 

(11.22) 

(11.23) 


where  A  is  the  total  area  of  the  domain  and  g(s)  is  some  function  of  the  distance 
along  the  domain  boundary.  Equation  11.21  results  from  the  Divergence  Theorem, 
which  states  that  the  volume  integral  of  the  divergence  of  a  vector  field  is  equal  to 
the  flux  of  the  vector  field  through  the  surface  of  the  volume.  The  two-dimensional 
analog  is  that  the  surface  integral  of  the  divergence  of  a  surface  vector  field  is  equal 
to  the  flux  of  the  vector  field  through  the  boundary.  In  mathematical  language, 


u  dA  =  u  ■  n  ds. 

J  dQ 


(11.24) 
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Combining  equation  11.24  with  II.  10  gives: 

[  [  V  •  u  dA  =  [  [  W2(p  dA  =  [  u  ■  n  ds  =  [  g(s )  ds.  (11.25) 

J  Jn  J  Jn  Jdci  Jdci 

The  left-hand  side  of  equation  11.21  represents  the  average  of  the  Laplacian  of  tpb 
over  the  domain.  Notice  that  the  interior  modes  of  ip  satisfy  equation  11.14,  while  the 
boundary  modes  of  y>  are  allowed  to  have  flow  through  the  boundary  (n  •  V<pb  ±  0). 
The  normal  flow  across  open  boundary  segments,  n  •  u\^n,  is  given  by  measured 
data,  if  available.  Other  data,  such  as  current  data  from  numerical  models,  may 
also  be  used  to  specify  normal  boundary  flow.  To  solve  equations  11.21  -  11.23,  the 
boundary  function  g(s)  is  expanded  in  terms  of  a  set  of  basis  functions.  Any  set 
of  basis  functions  can  be  used,  though  one  convenient  set  is  the  set  of  Fourier  basis 
functions  gi(s)  =  {cos  ^p,sin  ^p},  where,  again,  s  is  the  distance  along  the  domain 
boundary,  l  is  the  total  length  of  the  boundary,  and  i  is  a  basis  index.  The  boundary 
mode  process  is  similar  to  solving  a  nonhomogeneous  ODE,  whereby  the  solution 
to  the  homogeneous  problem  is  found  first,  and  then  a  particular  solution  to  the 
nonhomogeneous  problem  is  found.  In  the  case  of  OMA,  the  interior  modes  satisfy 
the  homogeneous  problem  (with  no  boundary  flow),  and  the  boundary  modes  satisfy 
the  nonhomogeneous  problem  (with  boundary  flow  specified  by  measured  data).  By 
the  principle  of  superposition,  the  sum  of  the  general  solution  and  the  particular 
solution  is  also  a  solution.  More  information  on  the  boundary  modes  is  available  in 
[7]  and  [8]. 

With  these  boundary  modes  added,  the  full  current  field  (including 
current  through  open  boundary  segments)  is  reconstructed  as: 

oo  oo  oo 

u  =  £  a*  (fe  x  VA)  +  £  atv.fi,  +  £  atv fit  (11.26) 

2=1  2=1  2=1 

Since  the  exact  solution  of  the  PDE  will  involve  an  infinite  number  of 
modes,  some  stopping  criteria  is  specified  in  the  mode  generation.  This  is  done  by 
specifying  a  mode  scale.  All  modes  on  a  scale  less  than  the  stopping  scale  are  thrown 
out.  More  detail  on  the  scale  and  stopping  criteria  is  contained  in  [7]. 
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2.  Unstructered  Grid 

In  contrast  to  the  total  currents  method,  the  OMA  method  is  not  calculated 
on  a  rectangular  grid.  The  PDE  is  instead  solved  on  an  unstructured  triangular  mesh. 
Rectangular  grids  can  be  problematic  for  solving  PDEs  on  arbitrary  non-rectangular 
domains.  If  greater  grid  resolution  is  desired  (along  a  coastline  for  example),  a  rect¬ 
angular  mesh  must  be  subdivided  across  the  entire  domain,  leading  to  vast  increases 
in  processing  time  and  required  resources.  Furthermore,  a  complex  coastline  approx¬ 
imated  by  rectangles  might  have  problems  if  a  no-flow  boundary  condition  is  applied 
on  the  artificial  ‘corners’  of  the  coast.  It  should  be  noted  that  in  our  case,  no  no-flow 
boundary  conditions  were  applied  in  the  rectangular  grid  (totals  method).  In  that 
case,  the  currents  were  generated  from  a  least-squares  fit  to  the  radial  measurements, 
and  any  current  vectors  falling  on  the  coastline  (from  erroneous  radial  measurements, 
for  example)  were  masked  out.  Currents  next  to  the  coastline,  however,  are  allowed 
to  have  current  vectors  that  appear  to  flow  through  the  coastline.  Figure  7  shows  an 
example  of  a  grid  mesh  and  an  unstructured  triangular  mesh  applied  near  a  coastline 
point. 

With  an  unstructured  triangular  mesh,  however,  higher  resolution  is  available 
in  the  areas  where  it  is  needed  without  unecessarily  increasing  resolution  across  the 
entire  domain.  In  addition,  no-flow  boundary  conditions  can  be  applied  in  a  more 
realistic  manner. 

3.  Fitting  Current  Data  to  OMA  Modes 

Once  the  modes  are  calculated,  any  current  field  (up  to  the  scale  of  modes 
that  were  kept)  can  be  represented  as  a  linear  combination  of  these  modes.  The 
linear  combination  is  written  in  equation  11.26.  The  coefficients  of  the  modes  in  this 
linear  combination  are  found  by  fitting  the  existing  current  measurements  to  the 
OMA  modes.  The  currents  are  fit  using  the  u  and  v  components  of  the  total  currents 
generated  from  the  radial  measurements  or  using  the  radial  measurements  directly. 
Using  radial  measurements  directly  to  find  the  OMA  mode  coefficients  is  generally 
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Longitude 


(a) 


(b) 


Figure  7.  Example  of  Coastline  Grid  around  Point  Pinos  at  the  souther  end  of  Mon¬ 
terey  Bay  using  (a)  Rectangular  Grid,  and  (6)  Unstructured  Triangular  Mesh.  Both 
panes  show  the  flow  conditions  on  the  domain  boundary.  Arrows  indicate  allowed 
direction  of  allowed  flow. 

preferable  because  it  avoids  the  additional  error  introduced  during  the  radials  to 
totals  calculation.  In  fitting  the  modes  to  existing  current  data,  the  mode  coefficients 
are  calculated  that  provide  the  minimal  error  between  modal  currents  and  currents 
measured  by  radar.  This  is  done  by  minimizing  a  cost  function.  The  cost  function 
is  defined  as  the  square  root  of  the  sum  of  the  squared  differences  between  modal 
currents  and  measured  currents,  although  other  cost  functions  could  be  used.  One 
adjustment  to  the  cost  function  used  in  the  processing  of  this  thesis  is  to  add  a 
penalty  for  large  mode  coefficients.  This  prevents  modal  currents  from  becoming 
unrealistically  large  in  areas  where  there  are  no  current  measurements  to  constrain 
them.  In  areas  where  measured  data  is  sparse,  modal  currents  can  become  quite  large 
without  affecting  the  difference  between  modal  currents  and  measured  currents.  The 
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large  coefficient  penalty  term  in  the  cost  function  reduces  the  tendency  of  the  fitted 
currents  to  become  unreasonably  large  in  areas  where  data  measurements  are  sparse. 
Additionally,  weights  may  be  introduced  to  place  more  importance  on  certain  data, 
although  this  is  not  done  in  this  thesis. 

The  general  cost  function  used  to  find  mode  coefficients  is: 


c  = 


M 

/ 

Wr 

ml 

\  m=  1 

N 
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^  '  f'r. 

\n=l 


U' 


M  * 

Z  n= 1 


(11.27) 


where  M  is  the  total  number  of  measurements,  N  is  the  total  number  of  modes,  un(xm ) 
is  the  current  at  location  m  associated  with  the  nth  mode,  r  is  the  unit  vector  from  the 
current  location  to  the  radar  site,  and  urm  is  the  measured  radial  current  at  location 
rn  in  the  direction  of  f.  To  minimize  this  function,  note  that  minimizing  the  square 
of  this  cost  function  is  equivalent  to  minimizing  the  cost  function  itself.  Taking  Aki 


0a„ 


and  setting  equal  to  zero  gives 


d(C2) 

da„ 
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=  E 

m= 1  L 
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m  >  (o^?7,/U77,(xm))  *  xm  ^mj 


+  M  Kjj  otj  —  0 


(11.28) 


for  all  coefficients  a3.  This  is  a  linear  set  of  M  equations  in  N  unknowns.  The  system 
can  be  solved  exactly  if  N  measurements  exist  (M  =  N ),  or  can  be  approximated 
with  a  least  squares  solution  if  more  than  N  measurements  exist  (M  >  N).  An 
overdetermined  system  is  preferable  to  an  exact  solution,  since  the  measurements 
contain  noise  which  can  presumably  be  filtered  out  if  many  more  measurements  than 
modes  exist  (M  S>  N). 

To  ensure  that  equation  11.28  is  indeed  a  minimum  of  the  cost  function  equation 
11.27,  the  second  derivative  is  examined. 

d(C2)2  M 


dal 


=  £  [2»£(« 


j  \xm) 


m—  1 


-\~  M  hij 


(11.29) 


Assuming  that  the  large  coefficient  penalty,  k,  and  the  weights,  W]n  are  both  posi¬ 
tive,  equation  11.29  is  everywhere  positive  and  equation  11.28  is  guaranteed  to  be  a 
minimum  of  the  cost  function  11.27. 
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E.  PARTICLE  TRAJECTORIES 


Once  the  current  field  is  known,  forward  particle  trajectories  can  be  calculated. 
A  particle  trajectory  describes  the  path  a  particle  takes  when  moving  with  a  velocity 
described  by  the  current  velocity  field.  The  trajectories  are  calculated  by  solving  the 
ordinary  differential  equation  (ODE): 

r!  nr* 

—  =  u(t,x(t)),  (11.30) 

where  x  is  the  position  vector  composed  of  latitude  and  longitude  coordinates  and 
u(t,x(t))  is  the  velocity  vector  composed  of  latitudinal  and  longitudinal  velocities. 

Many  techniques  are  available  for  numerical  solution  of  ODEs.  The  trajectory 
equation  11.30  is  solved  using  the  Runge-Kutta  fourth-order,  fifth-stage  method  dur¬ 
ing  this  investigation.  The  Runge-Kutta  fourth-order,  fifth-stage  method  for  solving 
ordinary  differential  equations  is  a  one-step,  multi-stage  method.  It  is  a  one-step 
method  in  that  to  find  the  value  of  x(ti+ 1),  only  data  from  the  previous  point  x{tj)  is 
used.  The  Runge-Kutta  fourth-order,  fifth-stage  method  uses  the  following  scheme: 
let  =  f(t,y)  be  an  ordinary  differential  equation  with  initial  condition  y(to)  =  yo . 
The  solution  to  this  initial  value  problem  can  be  approximated  numerically  by: 


Vn+i  —  Vn  +  h{b\k\  +  b2k2  +  b^ks  +  64A4  +  65^5),  (11.31) 

where 

k\  —  f(tn,  yn)  (11.32) 

k2  =  f (tn  +  c2h,  yn  +  h(a21ki ))  (11.33) 

^3  =  f(tn  +  c2>h-,  yn  +  h{a^\k\  +  ^32^2))  (11.34) 

k\  =  f(tn  +  C4I1,  yn  +  h(a4iki  +  a42k2  +  <3.43^3))  (11.35) 

^5  —  f  (tn  +  C5/1,  yn  +  /f(<7r,i/q  +  052^2  +  «53^3  +  ^54^4))-  (11.36) 


The  coefficients  at ,  b, ,  Ci  are  given  by  the  particular  fourth-order  Runge-Kutta 
algorithm  being  used  and  h  is  the  step  size.  The  fifth  stage  is  used  to  bound  the 
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error  of  the  approximation  yn+\.  The  error  can  be  controlled  by  examining  the  fifth 
stage  and  varying  the  step  size  h  to  keep  the  local  truncation  error  below  a  specified 
bound. 

In  the  case  of  this  thesis,  equation  11.30  is  the  ODE  that  is  solved  to  make 
particle  trajectories.  Let  xn  be  the  particle  position  at  time  tn.  Then  the  particle 
position  at  the  next  time  step  tn+ 1  is  approximated  by: 


xn+i  —  xn  +  h(b\k\  +  b2k2  +  +  64^)4  +  b3k§),  (11.37) 

where 

k1  =  u(tn,xn)  (11.38) 

k2  =  u(tn  +  c2h,  xn  +  h(a2ikl))  (11.39) 

k3  =  u(tn  +  c3h,  xn  +  h{a3\k\  +  a32k2))  (11.40) 

k4  =  u(tn  +  C4/1,  xn  +  h(a4iki  +  a42k2  +  ch43k3y)  (11.41) 

^5  =  u{tn  +  C5/1,  xn  +  h(a5iki  +  a32k2  +  a33k3  +  054^4)).  (11.42) 

tt(t,£c(t))  is  the  current  velocity  field. 


In  the  numerical  solution  of  many  ODEs,  evaluations  of  the  function  ^  = 
f(t,  y )  are  a  computational  expensive  step.  In  this  thesis,  evaluating  the  derivative  of 
the  position  vector  simply  involves  a  lookup  in  the  velocity  data  field,  which  contains 
velocity  information  at  specific  grid  points.  When  the  value  of  the  derivative  is 
needed  at  a  location  not  specified  in  the  current  fields  (any  location  other  than  the 
grid  points),  interpolation  is  done  to  approximate  the  currents  at  the  desired  points. 
If  the  current  field  is  the  output  of  the  totals  current  method,  then  the  field  is  defined 
on  a  regular  latitude/longitude  grid.  Simple  bilinear  interpolation  in  space  is  used  to 
interpolate  between  the  nearest  four  latitude/longitude  points  where  the  currents  are 
known.  If  the  current  field  is  the  output  of  the  OMA  method,  then  the  currents  are 
defined  on  a  triangular  mesh.  Currents  are  assumed  to  be  constant  across  a  triangle, 
and  if  a  current  measurement  is  needed  at  a  grid  point  that  is  not  the  center  of  a 
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triangle,  then  the  current  value  for  whatever  triangle  the  grid  point  falls  within  is 
returned.  This  is  nearest  neighbor  interpolation.  For  both  the  Totals  method  and 
the  OMA  method,  bilinear  interpolation  in  time  is  performed  last  to  find  current 
values  at  times  not  specified  in  the  current  field. 

Bilinear  interpolation  is  not  a  highly  accurate  interpolation  scheme.  Higher 
order  interpolation  methods  are  available  that  give  greater  accuracy  in  interpolating 
functions.  However,  the  sampling  error  inherent  in  the  current  fields  (error  on  the 
order  of  10%)  usually  will  dominate  any  error  resulting  from  interpolation,  making 
higher  order  interpolation  methods  unnecessary  [3]. 
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III. 


METHODS 


A.  MODEL  VALIDATION  AND  SST  COMPARISONS 

Figure  8  shows  a  flow  chart  of  the  processing  and  validation  procedure.  The 
SST  comparisons  formed  the  main  test  of  the  thesis  hypothesis.  SST  was  compared  in 
three  different  ways:  advected  SST  comparison,  static  SST  comparison,  and  random 
SST  comparison.  The  advected  method  of  SST  comparison  compared  SST  according 
to  the  hypothesized  model.  Static  and  random  comparison  methods  were  meant  as 
baselines.  In  each  case,  absolute  SST  and  SST  anomaly  were  used  in  the  comparison. 

Statistics  were  calculated  on  the  SST  differences  between  all  the  points  in  the 
comparison.  One  statistic  used  in  the  comparisons  is  the  root  mean  square  (RMS) 
difference.  The  RMS  is  defined  as: 


RMS* 


dif  ference 


\ 


J2(Tk{to,x0,y0) 
k= 1 


Tk(tuxi,yi)f 


n 


(HIl) 


That  is,  the  RMS  difference  is  the  square  root  of  the  average  of  the  squared  temper¬ 
ature  differences. 

In  the  advected  method,  SST  values  at  the  arrival  points  are  compared  with 
SST  values  at  the  departure  points.  Between  images,  many  points  are  compared.  An 
illustration  of  the  advected  comparison  is  shown  in  Figure  9. 

In  the  static  comparison  method,  SST  values  at  the  arrival  points  in  the  sec¬ 
ond  image  were  compared  against  SST  values  at  those  same  locations  in  the  first 
image.  Again,  statistics  were  calculated  on  the  SST  differences  of  all  the  points  in 
the  comparison.  An  illustration  of  the  static  comparison  method  is  shown  in  Figure 
10. 

In  the  random  comparison  method,  SST  values  at  the  arrival  points  in  the 
second  image  were  compared  against  random  points  within  the  same  image.  An 
illustration  of  the  random  comparison  method  is  shown  in  Figure  11. 
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Figure  8.  Flow  chart  of  SST  comparison  procedure. 
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Figure  9.  Example  illustration  of  comparison  of  advected  points  between  SST  images. 
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Figure  10.  Example  illustration  of  static  comparison  of  points  between  SST  images. 
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Figure  11.  Example  illustration  of  comparison  of  random  points  within  a  SST  images. 

In  testing  the  hypothesized  model,  all  possible  combinations  of  SST  images 
within  the  case  study  are  compared  with  each  other.  Each  comparison  method  (ad- 
vected,  static,  and  random)  was  used  in  the  comparison  between  each  pair  of  images. 
If  the  advected  method  is  a  good  predictor  of  SST,  then  the  RMS  statistic  of  the 
advected  comparisons  will  be  closer  to  zero  than  the  RMS  statistic  of  the  other  com¬ 
parison  methods. 

B.  DATA  USED 
1.  SST  Data 

As  mentioned  previously,  we  used  the  NOAA  high  resolution  AVHRR  SST 
dataset  in  this  thesis.  All  single-pass  POES  SST  images  were  downloaded  for  the  case 
study  periods  at  the  highest  resolution  available,  approximately  1.5km.  Data  from 
each  satellite  pass  within  the  case-study  periods  were  downloaded  from  the  NOAA 
Coast  Watch  Browser  in  MATLAB  format.  NOAA  Coast  Watch  data  is  available  at 
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http://coastwatch.pfel.noaa.gov/.  Each  SST  pass  was  analyzed  for  cloud  cover.  If  a 
pass  had  relatively  low  cloud  cover  within  the  CENC  domain,  it  was  flagged  to  be  used 
in  the  analysis.  These  good  images  were  used  because  of  the  high  probablility  that 
independent  SST  measurements  from  the  satellites  will  be  available  for  the  location 
we  need.  If  we  included  every  SST  image,  including  the  cloudy  ones,  then  few  SST 
measurements  would  be  available  for  comparison.  After  SST  images  were  downloaded 
and  good  images  selected  for  each  case  study,  the  images  were  masked  using  the  CENC 
area. 

2.  HF  Radar  Data 

As  mentioned  above,  the  HF  radar  data  used  in  this  thesis  originates  with  the 
California  Ocean  Currents  Monitoring  Program  (COCMP).  This  data  is  available  at 
http://www.cencalcurrents.org/.  Data  was  downloaded  as  individual  radial  hies  for 
every  hour  within  the  case  study  periods.  A  single  hie  contained  radar  data  for  one 
hour  from  a  single  radar  site. 

C.  CASE  STUDIES 

To  test  the  hypothesis  that  current  data  could  be  used  to  improve  SST  pre¬ 
dictions,  test  cases  were  needed. 

1.  CENC  Domain 

Our  area  of  study  is  the  coastal  ocean  of  central  California  from  Point  Sur  to 
Point  Arena  and  out  to  sea  approximately  100  nautical  miles.  This  area,  termed  the 
CENC  domain,  is  the  approximate  limit  to  the  HF  Radar’s  spatial  coverage.  High 
resolution  satellite  SST  data  (approximate  resolution  is  1.5  km)  is  available  for  this 
region  from  the  POES  AVHRR. 

2.  Dates  of  Study 

Two  time  periods  were  chosen  for  our  study.  These  case  study  periods  were 
chosen  because  of  relative  abundance  of  both  HF  Radar  data  and  clear  satellite  SST 
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images.  The  limiting  factor  in  case  study  selection  is  the  number  of  clear  SST  images 
available  for  comparison.  The  SST  data  from  these  case  study  periods  contain  a  high 
number  of  very  clear  images,  necessary  for  comparison  of  advected  SST  pixels  with 
SST  measured  from  satellites.  The  case  study  periods  were  October  17-28,  2006  and 
January  9-29,  2007. 

In  the  October  2006  case  study,  there  are  a  total  of  43  POES  satellite  passes 
with  an  average  cloud  cover  of  47.5%  within  the  CENC  domain.  Of  these  43  satellite 
passes,  11  passes  were  flagged  to  be  used  in  the  study.  These  11  images  have  an 
average  cloud  cover  of  8.9%.  In  the  January  2007  case  study,  there  are  a  total  of  74 
POES  satellite  passes  with  an  average  cloud  cover  of  57.6%  within  the  CENC  domain. 
Of  these  74  satellite  passes,  19  passes  were  flagged  to  be  used  in  the  study.  These  19 
images  have  an  average  cloud  cover  of  11.9%. 

Figure  12  shows  a  plot  of  the  normalized  radar  current  data  coverage  for  this 
region  for  the  time  periods  of  the  study.  The  colored  squares  in  Figure  12  represent  the 
percentage  of  time  a  radar  measurement  was  available  for  the  different  case  studies. 
The  white  outline  is  the  limits  of  the  CENC  region. 

3.  Stommel  Ocean  Model 

An  analytical  ocean  current  field  was  needed  to  compare  currents  generated  by 
the  Totals  and  OMA  method  to  an  analytical  solution.  While  any  vector  field  with 
no-flow  boundary  conditions  could  be  used,  the  Stommel  ocean  model  was  chosen 
because  of  its  familiar  use  in  oceanography.  The  Stommel  ocean  model  was  first  in¬ 
troduced  by  Henry  Stommel  to  study  the  causes  of  a  phenomenon  known  as  westward 
intensification  [9] .  Westward  intensification  refers  to  the  tendency  of  large-scale  ocean 
circulation  to  form  fast,  narrow,  and  deep  currents  on  the  western  edge  of  ocean  basins 
and  slow,  wide,  and  shallow  currents  on  the  eastern  edge  of  ocean  basins.  The  model 
assumes  a  rectangular  ocean  with  wind  stress,  Coriolis,  gravitational,  and  frictional 
forces.  The  component  of  flow  normal  to  any  of  the  ocean  boundaries  is  necessar¬ 
ily  zero.  Wind  stress  forces  are  defined  to  be  purely  in  the  east-west  direction  as, 
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(a)  (6) 

Figure  12.  Radar  data  coverage  for  the  central  California  region,  a.  October  17-28, 
2006  b.  January  9-29,  2007 

r  =  —  cos  westward  in  the  bottom  part  of  the  ocean,  and  eastward  in  the  top 
part  of  the  ocean.  With  these  assumptions,  the  model  can  be  represented  by  the 
equation: 


=  f)-  (-) 

JMs  a  stream  function  which  describes  the  current  velocity,  rG  is  the  maximum  am¬ 
plitude  of  wind  stress,  7  is  the  bottom  friction  coefficient,  f3  is  the  variation  of  the 
Coriolis  parameter  with  latitude,  and  L  is  the  length  and  height  of  the  domain.  A 
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steady-state  analytic  solution  to  this  equation  is  possible.  The  u  and  v  components 
of  the  steady-state  solution  are: 


u  = 


-^  =  (-l)(cieA‘I  +  C2^  +  C3)cos(T 


d 


v  =  —  =  (ci\ieAlX  +  c2A2e 


sm 


ny 


where 


(III.3) 

(III-4) 


Al,2  — 


C2 


7T  N 


£  ±t  (£)2  +  (y) 


(HI. 5) 
(ni-6) 

(HI. 7) 
(III.8) 


A  fictitious  rectangular  ocean  domain  was  created  and  used  to  simulate  steady 
state  Stommel-like  circulation.  To  provide  a  geophysical  context,  the  fictitious  ocean 
was  given  longitude  boundaries  of  —124°  West  to  —123°  West  and  latitude  boundaries 
of  36°  North  to  37°  North  (Figure  13).  It  should  be  noted  that  the  Stommel  ocean 
model  is  used  to  model  oceans  on  an  ocean  basin-wide  scale.  The  size  of  our  domain 
is  significantly  smaller  and  the  Stommel  model  does  not  have  a  physical  context  on 
these  scales.  The  Stommel  model  was  chosen  purely  for  its  analytic  current  field  and 
familiarity  to  those  in  the  field  of  oceanography. 

The  ocean  surface  flow  at  any  point  was  calculated  using  equations  III. 3  and 
III. 4.  Within  the  ocean  model,  a  subdomain  was  mapped  out  for  testing  of  the  HFR 
current  generation  methods.  This  subdomain  had  open  boundaries  on  all  sides  except 
the  southern  boundary,  which  was  shared  with  the  larger  ocean  domain  and  is  a  closed 
boundary. 
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Figure  13.  Stommel  Ocean  model  and  sample  HF  Radar  domain 

On  the  southern  boundary,  five  radar  sites  were  placed  at  regular  intervals. 
These  sites  were  used  to  sample  the  current  velocity  field  based  on  typical  resolutions 
of  HF  Radars.  For  each  radar  site,  a  radial  grid  of  points  was  created.  The  radial 
grid  points  were  spaced  at  5°  angles  between  radial  ‘spokes’  and  3  kilometers  distance 
between  radial  range  rings.  This  spacing  is  typical  of  the  radar  measurement  resolu¬ 
tion  seen  on  the  Central  California  coast.  At  these  radial  grid  points,  the  Stommel 
field  was  sampled  and  the  component  of  the  total  current  velocity  in  the  direction  of 
the  bearing  line  between  radar  site  and  grid  point  was  recorded.  The  component  of 
the  velocity  along  the  bearing  line  is: 

r  =  |w|  cos  (0  —  6*),  (III. 9) 


35 


where  |w|  is  the  magnitude  of  the  current  velocity  sampled  from  the  Stommel  field 
using  equations  III. 3  and  III. 4,  0  is  the  heading  from  the  radial  grid  point  to  the 
radar  site,  and  9  is  the  direction  of  the  current  vector.  A  positive  radial  component 
points  from  the  radial  grid  point  to  the  radar  site.  All  angles  are  measured  in  the 
traditional  mathematical  sense  of  counterclockwise  from  East. 

After  sampling  of  the  Stommel  current  field,  each  radar  site  had  measurements 
of  the  radial  current  components  at  the  locations  of  the  radial  grid  points.  A  random 
error  was  added  to  each  radial  measurement  to  simulate  instrument  or  sampling 
noise.  The  random  error  was  uniformly  distributed  from  -10%  to  +10%.  In  addition, 
a  certain  percentage  of  random  radial  measurements  were  thrown  out  to  represent 
missing  data.  The  error  and  missing  data  were  added  to  attempt  to  represent  physical 
reality,  which  contains  erroneous  and  missing  data.  These  radial  measurements  were 
then  processed  back  to  the  current  field  using  either  the  Totals  Method  or  the  OMA 
Method.  A  time-series  of  the  Stommel  current  field  was  built  representing  15  days 
of  data  with  a  sampling  every  hour,  for  a  total  of  360  measurements.  While  each 
sampling  of  the  Stommel  current  field  sampled  a  steady-state  field,  the  current  field 
resulting  from  the  Totals  or  OMA  method  was  different  for  each  hour’s  data,  due 
to  the  random  error  and  missing  measurements  introduced.  The  OMA  and  Totals 
method  were  analyzed  under  these  differing  conditions  to  see  which  method  proved 
more  robust  to  erroneous  or  missing  data. 

D.  MATLAB  PROCESSING 

In  processing  the  current  data,  we  used  the  ‘HFR  Progs’  MATLAB  toolbox 
developed  by  Dr.  David  Kaplan.  ‘HFR  Progs’  is  a  toolbox  for  the  processing,  view¬ 
ing,  and  analysis  of  HF  Radar  Currents.  The  latest  version  is  available  at 
https://cencalarchive.Org/~ cocmpmb/COCMP-wiki.  It  provides  functionality  to  pro¬ 
cess  radial  current  measurements  to  total  currents  and  particle  trajectories  using  the 
totals  method,  as  well  as  functionality  to  compute  OMA  modes  and  fit  radial  data 
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to  OMA  modes  and  calculate  trajectories  based  on  the  OMA  fits.  All  processing  was 
completed  using  either  MATLAB  7.0  or  MATLAB  7.4. 

1.  Processing  of  Radial  Data 

Radial  current  data  was  obtained  from  the  COCMP  Web  Archive  for  the  two 
case  studies.  Each  radar  data  file  contains  one  hour’s  worth  of  radial  surface  current 
for  one  radar  site.  Before  the  radial  files  are  processed  to  Total  currents  or  fit  to 
OMA  modes,  they  are  cleaned.  The  radial  measurements  are  cleaned  by  removing 
any  measurements  greater  than  100  cm/sec,  which  is  a  reasonable  upper  bound  on 
current  velocities  for  the  CENC  region  [10].  All  radial  measurements  from  each  radar 
site  are  also  cleaned  using  a  masking  polygon  unique  to  each  radar  site.  The  masking 
polygon  for  each  site  represents  the  area  of  reasonable  data  coverage.  Occasionally, 
the  radar  will  report  data  from  areas  which  are  unreasonable,  such  as  current  data 
lying  over  land.  If  the  radial  data  file  contains  measurements  that  he  outside  the 
masking  polygon,  these  measurements  will  be  removed  in  this  masking  step.  Finally, 
radial  data  was  interpolated  to  fill  in  gaps.  Interpolation  was  accomplished  using  the 
closest  values  in  the  bearing  direction  (filling  in  a  maximum  of  two  missing  bearing 
bins)  and  in  the  range  direction  (filling  in  a  maximum  of  one  missing  range  bin). 

2.  Processing  to  Totals  Currents 

Cleaned,  masked,  and  interpolated  radial  current  measurements  were  then 
processed  to  total  currents  using  the  least-squares  fit  described  in  Chapter  II,  section 
C.  Currents  were  generated  for  every  hour  of  each  day  of  the  case  study  time  periods. 
The  currents  were  generated  on  a  regular  latitude/longitude  grid  with  two  and  a 
half  kilometer  spacing,  which  is  comparable  to  the  native  resolution  of  the  radar 
observations  [3].  In  order  to  make  a  total  current  measurement  at  a  grid  point,  it  was 
specified  that  at  least  three  radial  measurements  from  at  least  two  different  radar  sites 
were  required.  All  radial  measurements  within  three  kilometers  of  a  grid  point  were 
used  in  each  total  current  generation.  With  this  search  radius,  it  is  possible  that  a 
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single  radial  measurement  may  have  been  used  in  making  total  current  measurements 
at  multiple  grid  points  (i.e. ,  neighboring  grid  results  are  not  completely  independent). 
After  currents  were  generated  on  the  grid  for  each  hour,  the  total  currents  were 
cleaned  and  masked.  The  currents  were  cleaned  by  removing  any  current  vector  with 
a  speed  above  100  cm/sec.  Total  current  vectors  were  cleaned  further  by  removing 
any  current  vector  with  a  GDOP  error  above  1.5,  where  GDOP  was  calculated  using 
the  two  most  orthogonal  radial  measurements  (error  flag  is  ‘GDOPMaxOrthog’  in  the 
HFR  Progs  MATLAB  Toolbox).  A  GDOP  error  less  than  1.5  corresponds  to  an  angle 
separation  of  between  71°  and  109°.  Any  current  measurements  remaining  after  the 
mask  was  applied  were  generated  from  radial  measurements  with  angle  separations 
in  that  range. 

After  cleaning,  the  total  currents  were  again  masked  with  a  polygon  outlining 
the  CENC  domain  to  ensure  that  no  outlying  current  vectors  remained.  These  cur¬ 
rents  were  then  interpolated  in  space  and  time  to  fill  in  any  holes  in  the  space/time 
vector  current  grid.  The  interpolation  method  uses  bilinear  interpolation  in  space 
and  time,  and  then  takes  the  average  of  the  two. 

3.  OMA  Modes  Processing 

The  HFR  Progs  Toolbox  was  used  to  generate  OMA  modes  on  the  CENC  and 
the  Stommel  domain.  This  toolbox  makes  heavy  use  of  MATLAB’s  PDE  Toolbox  to 
generate  the  adaptive  mesh  and  to  solve  the  PDE  eigenvalue  problems  of  equations 
11.13  and  11.14.  MATLAB’s  PDE  Toolbox  uses  a  finite  element  method  to  numerically 
solve  eigenvalue  PDEs  such  as  this.  Certain  parameters  are  required  in  the  modes 
calculation.  The  modes  were  generated  with  a  cutoff  scale  of  five  kilometers.  A 
seed  value  for  the  estimated  number  of  triangles  in  the  unstructured  mesh  was  given 
as  10,000  triangles.  The  method  first  generated  a  triangular  mesh  to  use  in  the 
mode  generation.  The  MATLAB  toolbox  uses  a  Delaunay  triangulation  algorithm  to 
generate  the  mesh.  Figures  14  and  15  show  the  triangular  mesh  generated  for  the 
entire  domain  and  for  a  closeup  of  the  Monterey  Bay. 


38 


Figure  14.  Triangular  Mesh  for  CENC  region. 
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Figure  15.  Close-up  of  Triangular  Mesh  for  the  Monterey  Bay. 
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Modes  were  then  generated  on  this  triangular  mesh.  This  resulted  in  84  bound¬ 
ary  modes,  576  Dirichlet  (divergence-free)  modes,  and  653  Neumann  (vorticity-free) 
modes  for  a  total  of  f,313  modes.  This  is  a  larger  number  of  modes  for  a  larger  area 
than  previously  attempted  with  HF  Radar  data.  The  largest  scale  modes  are  plotted 
in  Figures  16,  17,  and  18. 

Modes  were  also  calculated  for  the  hypothetical  Stommel  model  using  the  same 
parameters. 

4.  Fitting  Radial  Data  to  OMA  Modes 

For  fitting  radial  data  to  OMA  modes,  the  radial  files  were  prepared  in  the 
same  method  as  preparing  for  total  current  generation,  except  no  radial  interpolation 
was  performed.  A  k  value  of  10-3  was  used  for  the  coefficient  penalty.  The  reasoning 
for  this  k  value  is  discussed  in  chapter  IV,  section  A. 

5.  Particle  Trajectories 

The  particle  trajectories  were  calculated  in  a  similar  way  for  both  the  Totals 
and  OMA  methods.  Once  currents  were  available,  both  methods  used  MATLAB’s 
‘ode45’  function  to  solve  the  advection  differential  equation  11.30.  The  MATLAB 
‘ode45’  function  is  a  Runge-Kutta  (4,5)  ODE  solver  which  uses  Runge-Kutta  coeffi¬ 
cients  developed  in  [11], 
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(a) 


(b) 


(c) 


(d) 


Figure  16.  Dirichlet,  or  Divergence-free  modes,  a,  b,  c,  d  show  the  first,  second,  third, 
and  fourth  divergence- free  modes:  xpi,  ip2,  ^3,  t/h-  Also  shown  are  the  corresponding 
contribution  to  the  current  field  from  these  modes  in  black. 
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Figure  17.  Neumann,  or  vorticity-free  modes,  a,  b,  c,  d  show  the  first,  second,  third, 
and  fourth  vorticity-free  modes:  </?i,  </?2,  <^3,  </?4-  Also  shown  are  the  corresponding 
contribution  to  the  current  field  from  these  modes  in  black. 
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Figure  18.  Boundary  modes,  a,  b,  c,  d  show  the  first,  second,  third,  and  fourth 
boundary  modes:  tpb2,  ip\.  Also  shown  are  the  corresponding  contribution  to 

the  current  field  from  these  modes  in  black. 
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IV. 


RESULTS 


A.  STOMMEL  MODEL  CURRENTS 
1.  Current  Statistics 

Using  an  analytic  solution  and  calculating  Totals  and  OMA  currents  from 
this  solution  allowed  for  direct  comparison  of  the  currents  and  the  analytic  solution. 
Radial  measurements  were  sampled  from  the  (steady-state)  analytic  solution  and  used 
to  build  up  currents  using  the  Totals  and  OMA  method.  Current  magnitudes  for  the 
analytic  solution  in  the  Stommel  subdomain  range  from  0  cm/sec  to  80  cm/sec,  with 
typical  values  in  the  20  -  30  cm/sec  range.  Error  and  sampling  holes  were  introduced 
as  described  in  Chapter  III,  section  3.  The  percentage  holes  varied  from  0%  to  90% 
in  increments  of  10%.  The  resulting  current  fields  were  then  compared  against  the 
analytic  solution.  Statistics  were  calculated  for  the  magnitude  of  the  difference  vector 
of  the  Total  and  OMA  current  field  and  the  analytic  solution.  These  statistics  are 
presented  in  Tables  I  and  II.  The  RMS  statistic  is  the  RMS  of  the  magnitude  of 
the  difference  vector,  and  mean  error  is  calculated  as  the  mean  of  the  ratio  of  the 
magnitude  of  the  difference  vector  and  the  magnitude  of  the  analytic  velocity  vector. 


%  Radial  Holes 

%  Missing  Totals  Points 

Max  Difference  Vector 

RMS  Difference 

Mean  Error 

0% 

0.0% 

1.7  cm/s 

0.3  cm/s 

5.9% 

10% 

0.3% 

2.2  cm/s 

0.3  cm/s 

6.3% 

20% 

0.8% 

2.6  cm/s 

0.4  cm/s 

6.8% 

30% 

1.7% 

3.0  cm/s 

0.4  cm/s 

7.4% 

40% 

2.9% 

3.4  cm/s 

0.5  cm/s 

7.8% 

50% 

6.1% 

4.1  cm/s 

0.5  cm/s 

8.3% 

60% 

13.3% 

4.7  cm/s 

0.6  cm/s 

8.1% 

70% 

23.2% 

4.6  cm/s 

0.6  cm/s 

7.3% 

80% 

44.6% 

4.4  cm/s 

0.6  cm/s 

5.0% 

90% 

78.3% 

6.7  cm/s 

1.0  cm/s 

1.9% 

Table  I.  Statistics  for  Totals  Currents  calculated  on  Stommel  domain. 
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%  Radial  Holes 

%  Missing  Totals  Points 

Max  Difference  Vector 

RMS  Difference 

Mean  Error 

0% 

0% 

41.3  cm/s 

0.8  cm/s 

11.6% 

10% 

0% 

41.0  cm/s 

0.8  cm/s 

11.9% 

20% 

0% 

40.3  cm/s 

0.8  cm/s 

12.3% 

30% 

0% 

40.0  cm/s 

0.9  cm/s 

13.0% 

40% 

0% 

39.5  cm/s 

0.9  cm/s 

13.8% 

50% 

0% 

39.7  cm/s 

1.0  cm/s 

14.9% 

60% 

0% 

41.4  cm/s 

1.1  cm/s 

17.0% 

70% 

0% 

47.0  cm/s 

1.2  cm/s 

19.4% 

80% 

0% 

60.4  cm/s 

1.5  cm/s 

25.5% 

90% 

0% 

134.0  cm/s 

2.6  cm/s 

44.9% 

Table  II.  Statistics  for  OMA  Currents  calculated  on  Stommel  domain. 


2.  Totals  Currents  Reconstruction 

Tables  I  and  II  reinforce  that  the  Totals  method  is  a  local  fit  method,  while  the 
OMA  method  is  a  global  fit  method.  It  is  seen  from  Table  I  that  the  Totals  currents 
are  a  good  approximation  for  the  analytic  currents  even  when  a  large  percentage  of 
the  radar  measurements  that  go  into  making  Totals  currents  are  removed.  This  is 
reflected  in  the  low  RMS  difference  and  mean  error  entries.  An  interesting  feature 
of  the  Totals  method  is  reflected  in  the  mean  error  column.  As  the  number  of  holes 
increases,  the  mean  error  of  the  measurements  also  increases,  but  then  starts  to 
decrease  as  the  percentage  of  radial  holes  grows  past  50%.  Due  to  the  spacing  and 
location  of  the  radar  sites,  the  final  current  measurements  become  concentrated  in 
the  lower  half  of  the  domain.  Even  when  90%  of  radial  measurements  are  removed, 
more  than  enough  radial  measurements  exist  to  generate  currents  in  the  southern 
portion  of  the  domain.  This  leads  to  a  very  good  fit  in  the  lower  part  of  the  domain, 
and  no  measurements  at  all  in  the  rest  of  the  domain,  hence  the  lower  mean  error. 
The  biggest  drawback  to  this  method  is  the  large  amounts  of  missing  data  that  result 
when  sparse  radial  coverage  is  available.  Figure  19  shows  two  reconstructed  current 
fields  using  the  Totals  method  on  the  Stommel  field  with  0%  and  90%  holes.  These 
are  single  time  slices  of  the  fields  which  give  rise  to  the  statistics  in  Table  I,  rows  one 
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and  ten.  The  reader  is  also  referred  back  to  Figure  13  in  Chapter  III,  section  C  for  a 
plot  of  the  analytic  Stommel  currents. 
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Figure  19.  Reconstruction  of  Stommel  current  held  using  Totals  method,  (a)  Current 
held  sampled  with  10%  error  and  0%holes.  (6)  Current  held  sampled  with  10%  error 
and  90%  holes. 


3.  OMA  Currents  Reconstruction 

In  contrast  to  the  Totals  method,  the  OMA  method  is  a  global  method.  All 
available  radial  measurements  are  used  in  the  cost  function  to  hnd  the  mode  coeffi¬ 
cients,  a:*,  that  best  ht  the  modes  to  available  measurements.  Because  the  modes  are 
dehned  over  the  entire  domain,  reducing  the  number  of  radial  measurements  will  not 
reduce  the  coverage  of  the  OMA  currents.  As  can  be  seen  in  the  OMA  statistics  table, 
however,  reducing  the  number  of  radial  measurements  does  reduce  the  accuracy  of  the 
OMA  ht.  These  large  differences  are  present  in  areas  with  few  radial  measurements. 
Even  for  relatively  dense  and  uniform  radial  coverage,  the  OMA  method  can  lead  to 
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large  differences  for  a  few  current  vectors,  particularly  along  the  edge  of  the  domain. 
Figure  20  shows  a  reconstructed  OMA  current  held  for  the  same  Stommel  held  with 
0%  and  90%  holes.  These  are  the  single  time  slices  of  the  helds  which  give  rise  to 
the  statistics  in  Table  II,  rows  one  and  ten.  It  is  seen  that  the  domain  boundaries, 
particularly  the  southwest  corner,  are  hard  for  the  OMA  method  to  ht. 
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Figure  20.  Reconstruction  of  Stommel  current  held  using  OMA  method,  (a)  Current 
held  sampled  with  10%  error  and  0%holes.  (6)  Current  held  sampled  with  10%  error 
and  90%  holes. 


Figure  21  shows  the  contribution  to  the  OMA  currents  from  each  of  the  differ¬ 
ent  type  of  modes:  Neumann  (99),  Dirichlet  (-0),  and  Boundary  (<%).  The  full  current 
held  is  a  sum  of  these  three  parts.  The  modes  are  ht  to  the  held  with  0%  radal 
holes  and  10%  error  (the  hrst  row  in  Table  II).  From  these  hgures,  it  is  seen  that 
the  boundary  modes  are  contributing  most  to  the  errors  seen  in  the  lower  left  and 
lower  right  portions  of  the  Stommel  domain.  The  Dirichlet  and  Neumman  modes  do 
not  have  much  difficulty  in  these  corners  because  the  boundary  conditions  that  are 
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imposed  on  these  modes  mean  that  there  is  no  difference  between  open  and  closed 
boundary  segments.  The  Boundary  modes,  however,  do  “see”  a  difference  in  open  and 
closed  boundaries.  These  corners,  where  a  closed  boundary  meet  an  open  boundary 
can  be  difficult  to  fit  for  boundary  modes.  The  biggest  errors  in  the  OMA  fits  occur 
in  these  corners,  where  the  boundary  modes  have  difficulty  adjusting  for  the  change 
from  closed  to  open  boundaries. 

4.  Mode  Coefficient  Penalty 

A  study  was  done  on  the  k  parameter  for  the  Stommel  domain.  Small  k 
values  translate  into  small  penalties  for  large  mode  coefficients.  If  there  are  few 
holes  and  abundant  radial  coverage,  then  the  measurements  naturally  limit  the  mode 
coefficients  from  becoming  large  and  the  n  term  in  the  cost  function  is  unnecessary. 
If  there  are  large  holes  in  the  radial  data  coverage,  then  the  k  term  is  necessary  in  the 
cost  function  to  keep  mode  coefficients  and  currents  from  becoming  large  in  areas  of 
little  data  coverage.  However,  including  a  k  term  which  is  too  large  might  suppress 
mode  coefficients  below  their  accurate  value.  Of  course,  the  extent  of  data  coverage 
is  not  known  before  the  time  of  measurement,  so  a  value  for  k  should  be  chosen 
that  both  suppresses  unrealistically  large  mode  coefficients  while  allowing  those  mode 
coefficients  to  accurately  describe  available  measurements.  Figure  22  shows  the  RMS 
difference  between  OMA  currents  and  analytic  currents  plotted  against  the  logarithm 
of  k  for  varying  amounts  of  missing  data.  A  k  value  of  10-3  —  10-2  minimizes  the 
RMS  difference  of  fitted  currents  with  actual  analytic  currents  for  a  wide  range  of 
missing  data.  A  value  of  10-3  was  chosen  for  this  investigation.  These  results  agree 
with  [7],  where  an  analysis  of  the  k  parameter  was  conducted  by  comparing  OMA 
and  Totals  currents  for  cases  near  Bodega  Bay. 
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Figure  21.  Contributions  to  the  Stommel  currents  from  the  different  types  of  modes. 
a.  Dirichlet  (curl  or  vorticity)  modal  currents,  b.  Neumann  (divergence)  modal  cur¬ 
rents.  c.  Boundary  modal  currents. 
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Figure  22.  RMS  difference  of  OMA  currents  and  the  analytic  Stommel  solution  for 
varying  amounts  of  missing  data. 

B.  SST  COMPARISONS  AND  MODEL  TESTING 

Currents  and  trajectories  were  calculated  for  the  CENC  domain  for  the  Oc¬ 
tober  and  January  case  studies.  SST  comparisons  were  carried  out  between  all  high- 
quality  satellite  image  pairs  in  the  case  studies  (both  for  absolute  temperature  and 
temperature  anomaly)  as  described  in  Chapter  III,  section  A.  For  the  October  case 
study,  this  resulted  in  a  total  of  55  image  pair  comparisons.  For  the  January  case 
study,  this  resulted  in  171  image  pair  comparisons.  The  result  of  these  comparisons 
were  a  set  of  statistics  on  each  pair  of  SST  images  that  were  compared.  Figures  23 
to  30  show  the  RMS  of  the  difference  of  SST  points  for  the  difference  comparison 
schemes  against  the  difference  in  time  between  compared  SST  images. 

Figures  23  through  26  show  the  absolute  SST  comparisons  using  both  Totals 
and  OMA  currents  for  the  two  different  case  studies.  Figures  27  through  30  show  the 
SST  anomaly  comparisons  using  both  Totals  and  OMA  currents  for  the  two  different 
case  studies. 
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Figure  23.  RMS  of  the  difference  of  SST  points  using  advected,  static,  and  random 
comparisons.  Absolute  SST  is  compared,  and  particle  trajectories  are  calculated  using 
Totals  currents  for  the  January,  2007  case  study. 
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Figure  24.  RMS  of  the  difference  of  SST  points  using  advected,  static,  and  random 
comparisons.  Absolute  SST  is  compared,  and  particle  trajectories  are  calculated  using 
OMA  currents  for  the  January,  2007  case  study. 
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Figure  25.  RMS  of  the  difference  of  SST  points  using  advected,  static,  and  random 
comparisons.  Absolute  SST  is  compared,  and  particle  trajectories  are  calculated  using 
Totals  currents  for  the  October,  2006  case  study. 
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Figure  26.  RMS  of  the  difference  of  SST  points  using  advected,  static,  and  random 
comparisons.  Absolute  SST  is  compared,  and  particle  trajectories  are  calculated  using 
OMA  currents  for  the  October,  2006  case  study. 
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Figure  27.  RMS  of  the  difference  of  SST  points  using  advected,  static,  and  random 
comparisons.  SST  anomaly  is  compared,  and  particle  trajectories  are  calculated  using 
Totals  currents  for  the  January,  2007  case  study. 
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Figure  28.  RMS  of  the  difference  of  SST  points  using  advected,  static,  and  random 
comparisons.  SST  anomaly  is  compared,  and  particle  trajectories  are  calculated  using 
OMA  currents  for  the  January,  2007  case  study. 
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Figure  29.  RMS  of  the  difference  of  SST  points  using  advected,  static,  and  random 
comparisons.  SST  anomaly  is  compared,  and  particle  trajectories  are  calculated  using 
Totals  currents  for  the  October,  2006  case  study. 
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Figure  30.  RMS  of  the  difference  of  SST  points  using  advected,  static,  and  random 
comparisons.  SST  anomaly  is  compared,  and  particle  trajectories  are  calculated  using 
OMA  currents  for  the  October,  2006  case  study. 
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Interpretation  of  these  statistics  is  difficult,  and  definite  conclusions  on  a  ‘best 
method’  impossible.  In  these  SST  comparisons,  the  random  comparison  is  used  as 
a  ‘baseline’  comparison.  In  the  random  comparison,  SST  (or  SST  anomaly)  at  the 
arrival  points  of  the  trajectories  are  compared  against  another  random  SST  pixel 
within  the  image.  This  serves  as  a  measure  of  the  inherent  variability  within  the 
image.  In  the  absolute  comparisons  in  the  January  case  study  (Figures  23  and  24), 
the  static  and  advected  methods  show  a  general  upward  trend  as  the  time  separation 
between  SST  images  increases  for  both  OMA  and  Totals  trajectories.  There  does 
not  appear  to  be  a  clear  distinction  between  the  advected  and  static  methods  for 
the  absolute  temperature  comparison  in  January  until  the  time  difference  reaches 
approximately  8-10  days,  at  which  point  it  appears  that  the  static  method  more 
regularly  leads  to  smaller  RMS  values  than  the  advected  method,  although  the  points 
are  still  grouped  together.  Static  and  advected  RMS  values  remain  lower  than  random 
comparisons  until  the  time  difference  is  approximately  6-8  days,  at  which  point  all 
comparison  methods  become  intermingled. 

For  absolute  temperature  comparisons  in  October,  advected  and  static  meth¬ 
ods  show  similar  results,  starting  lower  than  the  random  comparisons,  rising  up  higher 
than  the  random  points  around  time  differences  of  3-5  days,  then  falling  again  and 
starting  to  rise.  The  reason  for  the  peak  at  3-5  days  is  unclear.  In  general,  it  appears 
that  the  static  method  gives  lower  RMS  values  than  the  advected  method,  although 
the  points  are  usually  grouped  close  together. 

In  the  SST  anomaly  comparisons,  the  approximations  of  the  advected  and 
static  methods  become  more  accurate.  Anomalies  were  compared  as  an  attempt  to 
remove  the  effects  of  large  scale  warming  or  cooling  of  the  SST  in  our  domain.  If  the 
entire  area  was  warmed  or  cooled  at  a  similar  rate  over  the  time  difference  between 
starting  and  ending  images,  then  it  was  hoped  that  taking  the  anomaly  would  account 
for  most  of  this  change.  In  the  January  time  period,  both  advected  and  static  methods 
appear  to  perform  (predict  SST)  better  than  a  random  sampling.  For  time  differences 
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up  to  3-4  days,  the  two  methods  are  indistinguishable  for  both  OMA  and  Totals 
trajectories.  After  3-4  days,  the  static  method  generates  lower  RMS  differences  than 
the  advected  method.  Both  methods  seem  to  do  better  than  the  random  comparison 
method.  The  static  method  seems  to  outperform  the  advected  method  for  all  values  of 
time  separation,  although  they  are  grouped  fairly  close  together.  In  the  October  time 
period,  the  static  and  advected  methods  show  a  similar  pattern  as  in  the  absolute 
comparisons  (rising  RMS  values  reaching  a  peak  around  4-5  days,  then  falling  and 
rising  again).  In  this  case,  the  static  and  advected  methods  stay  below  the  random 
comparison.  There  is  no  discernable  difference  between  static  and  advected  methods 
in  this  time  period. 

1.  Comparison  of  OMA  and  Totals  Advected  Methods 

If  the  advected  method  offered  improvements  over  the  other  comparison  meth¬ 
ods,  then  it  is  expected  that  the  present  SST  field  could  be  predicted  by  a  past  SST 
field  when  advection  of  surface  water  is  taken  into  account.  This,  of  course,  will 
only  work  when  accurate  surface  current  and  particle  trajectories  are  available.  As  a 
measure  of  performance  of  the  two  different  surface  current  measurements,  OMA  and 
Totals  Currents,  the  advected  SST  comparisons  using  both  OMA  and  Totals  trajec¬ 
tories  are  plotted  together  for  the  two  case  studies  in  Figures  31  to  34.  These  figures 
are  the  same  statistics  for  the  advected  method  presented  in  Figures  23  to  30.  They 
are  extracted  and  presented  on  the  same  figures  to  ease  comparison.  Both  OMA  and 
Totals  currents  give  rise  to  similar  performance  of  the  advected  method  for  all  case 
studies. 
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Figure  31.  RMS  of  the  difference  of  SST  points  using  advected  comparisons  for  both 
OMA  and  Totals  trajectories.  Absolute  SST  is  compared  for  the  January,  2007  case 
study. 
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Figure  32.  RMS  of  the  difference  of  SST  points  using  advected  comparisons  for  both 
OMA  and  Totals  trajectories.  SST  anomaly  is  compared  for  the  January,  2007  case 
study. 
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Figure  33.  RMS  of  the  difference  of  SST  points  using  advected  comparisons  for  both 
OMA  and  Totals  trajectories.  Absolute  SST  is  compared  for  the  October,  2006  case 
study. 


1.2 

1 


•  Advect  OMA 
■  Advect  Totals 


°  0.8- 
O) 

CD 

CO 
DC 


0.6 


0.4 


0.2 


t 

« 

* 


•: 


.i 


4  5  6 

Time  Difference  (days) 


9  10 


Figure  34.  RMS  of  the  difference  of  SST  points  using  advected  comparisons  for  both 
OMA  and  Totals  trajectories.  SST  anomaly  is  compared  for  the  October,  2006  case 
study. 
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2.  Calculations  on  the  Distance  of  Advected  Particles 

In  the  advected  comparison  method,  calculations  were  also  completed  on  the 
distance  individual  particles  were  advected.  Individual  particle  displacement  dis¬ 
tances  were  combined  into  a  RMS  of  advected  displacement  statistic.  These  statistics 
are  plotted  in  Figure  35.  Predictably,  as  time  increases,  particles  are  advected  further 
from  their  starting  position,  on  average. 

C.  RECONSTRUCTION  OF  SST  FIELD 

If  the  advected  model  is  a  good  prediction  of  future  SST,  then  some  method 
is  needed  to  reconstruct  the  full  SST  field  from  the  (relatively)  few  points  that  are 
advected  from  one  image  to  another.  Numerous  interpolations  methods  exist  for  scat¬ 
tered  data  interpolation.  Figure  36  shows  the  reconstructed  and  actual  temperature 
field  for  October  25,  2006,  21:36  GMT.  The  reconstructed  field  is  generated  using  the 
SST  points  advected  (with  Totals  currents)  from  October  17,  2006,  21:18  GMT  to 
October  25,  2006,  21:36  GMT.  The  reconstruction  is  accomplished  using  MATLAB’s 
‘griddata’  function,  which  utilizes  triangle  based  linear  interpolation. 

The  reconstructed  field  and  the  actual  SST  field  are  significantly  different.  Ma¬ 
jor  features  of  the  actual  field  (Figure  36d),  such  as  the  cold-water  plume  off  Pescadero 
and  Half  Moon  Bay  in  the  center  of  the  image  are  absent  in  the  reconstructed  image. 
The  reconstructed  image  might  be  showing  some  promise  in  that  locations  of  SST 
features  such  as  fronts  (areas  of  high  gradient)  appear  to  be  in  the  correct  locations, 
but  overall  the  reconstruction  is  of  low  quality.  A  better  approximation  would  be 
Figure  36c,  which  is  just  the  past  SST  field  with  no  advection  information. 

An  alternative  method  of  reconstructing  the  SST  field  is  to  identify  points 
where  the  SST  is  needed,  and  the  find  the  backwards  trajectory  from  that  point  to  a 
point  in  a  previous  image  where  SST  is  available.  This  can  be  done  using  the  same 
current  information  necessary  to  compute  forward  trajectories. 
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Figure  35.  Displacement  of  advected  particles.  Particles  from  both  case  study  periods 
are  shown  in  each  plot,  (a)  Displacement  of  particles  using  OMA  currents.  (6) 
Displacement  of  particles  using  Totals  currents. 
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Figure  36.  (a)  Predicted  SST  field  for  21:36  GMT,  October  25,  2006,  based  on  par¬ 
ticles  advected  from  21:18  GMT,  October  17,  2006.  (6)  Reconstructed  temperature 
field  based  on  advected  points  in  panel  (a),  (c)  Actual  SST  field  for  21:18  GMT, 
October  17,  2006  as  measured  by  satellite.  ( d )  Actual  SST  field  for  21:36  GMT, 
October  25,  2006,  as  measured  by  NOAA  satellite. 
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D.  RECOMMENDATIONS  FOR  FURTHER  STUDY 

1.  Refinement  of  Advection  Model 

It  is  seen  that  the  advection  SST  model  proposed  is  not  an  adequate  predictor 
of  the  present  or  future  SST  held.  It  appears  that  the  static  model,  where  the  SST  at 
a  point  is  simply  modeled  by  a  SST  measurement  at  that  same  point  in  the  past,  does 
just  as  well  at  predicting  a  future  SST  value.  The  static  model  does  not  utilize  any 
information  about  the  surface  current  held.  The  advected  model,  in  contrast,  contains 
no  information  on  the  past  SST  at  that  location.  Perhaps  a  model  which  includes 
information  on  both  the  advected  SST  and  static  SST  would  perform  well.  This  is  a 
regression  model  question.  What  is  the  ideal  (if  any)  combination  of  advected  and 
static  SST  to  be  used  in  predicting  future  SST? 

2.  Frontal  Advection 

Another  proposed  advection  model  concentrates  the  predictions  in  areas  of 
interest.  It  is  often  the  case  that  researchers  are  more  interested  in  certain  features 
of  the  SST  held  than  others.  These  features  are  typically  areas  of  high  temperature 
gradients,  called  fronts.  Instead  of  predicting  a  SST  held  based  on  the  advected 
locations  of  a  grid  of  points  spread  over  the  entire  domain,  the  held  could  be  predicted 
by  the  advected  locations  of  a  set  of  points  that  dehne  the  feature  of  interest.  The 
advected  locations  of  the  features  can  be  used  as  control  points  and  the  entire  image 
registered  based  on  these  control  points. 

3.  Lagrangian  Coherent  Structures 

The  method  of  calculating  particle  trajectories  by  integrating  the  surface  ve¬ 
locity  can  lead  to  signihcant  errors  in  the  trajectories  for  small  changes  in  the  velocity 
held,  such  as  those  added  by  measurement  error.  Studies  are  underway  which  seek  to 
calculate  the  locations  of  more  robust  structures  in  the  surface  velocity  held  which  are 
less  sensitive  to  errors  in  individual  current  measurements.  Some  of  these  structures, 
called  Lagrangian  Coherent  Structures  (LCS),  form  natural  barriers  to  huid  transport 
and  act  as  separators  of  areas  of  different  dynamics.  These  LCS  can  be  calculated 
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from  the  Finite-Time  Lyapunov  Exponent  (FTLE),  discussed  in  [12].  Figure  37  is  an 
example  of  FTLE  and  LCS  caculations  from  HF  radar  current  data  in  the  Monterey 
Bay.  Note  the  LCS  running  north-south  across  the  month  of  Monterey  Bay  divides 
the  flow  within  the  bay  from  the  flow  which  remains  outside  the  bay.  LCS  can  be 
utilized  in  SST  prediction  by  helping  to  guide  calculated  trajectories,  reducing  the 
impact  of  measurement  errors  on  the  overall  particle  trajectories. 

4.  Adaptive  Interpolation  Functions  for  Scattered  Data 
Interpolation 

A  different  proposed  model  is  related  to  scattered  data  interpolation.  Instead 
of  modeling  a  present  (but  incomplete)  SST  field  by  a  past  field,  the  present  incom¬ 
plete  field  might  be  filled  in  by  interpolation  of  the  measurements  that  exist  in  the 
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Figure  37.  FTLE  field  computed  from  HF  radar  velocity  data.  Curves  of  high  FTLE 
represent  time- varying  LCS.  Also  shown  is  HF  radar  velocity  field  at  the  given  time. 
From  [12]. 
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present  field.  Existing  interpolation  methods  might  interpolate  a  value  based  on  a 
weighting  function  applied  to  surrounding  measurements.  One  weighting  function 
is  the  Gaussian  weighting  function,  which  weighs  surrounding  measurements  equally 
based  on  the  negative  exponential  of  the  squared  the  distance  to  the  measurement. 
The  Gaussian  function  is  shown  in  Figure  38. 

Continuing  on  the  logic  that  the  surface  current  field  contains  information 
which  is  relevant  to  the  surface  SST  field,  the  interpolation  weighting  function  might 
be  modified  based  on  existing  current  measurements.  For  example,  the  weighting 
function  can  be  adaptively  stretched  to  weigh  measurements  which  lie  along  local 
current  vectors  more  than  measurements  which  he  perpendicular  to  the  local  current 
vectors.  This  interpolation  function  can  be  changed  for  each  location  on  the  map 
based  on  local  current  information.  Figure  39  shows  a  interpolation  weighting  function 
based  on  the  Gaussian  function  which  has  been  compacted  in  the  x  direction. 


65 


Figure  39.  Gaussian  weighting  function,  slimmed  in  the  x  direction. 

5.  Improvements  to  Current  Fits  from  OMA  Method 

Additional  work  remains  on  the  OMA  method.  Improvements  can  be  made 
to  fitting  OMA  modes  to  existing  current  measurements.  In  this  investigation,  the 
only  measurements  used  to  fit  OMA  modes  originate  from  the  CODAR  HF  radars. 
Other  current  measurements  exist  which  can  be  used  to  fit  OMA  modes,  although 
care  should  be  taken  to  ensure  that  the  ocean  depth  at  which  currents  are  measured 
are  compatible  across  different  measurement  instruments.  For  example,  current  in¬ 
struments  on  moored  ocean  buoys  can  be  used  in  OMA  fits,  as  well  as  currents  from 
tidal  gauges.  Currents  near  land  are  hard  for  the  HF  radars  to  measure,  but  can  play 
a  large  part  in  current  circulation  patterns,  especially  in  bays  and  other  inlets. 
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APPENDIX 


Data  from  the  satellite  passes  in  Tables  III  and  IV  were  used  in  this  thesis. 


Date 

Time(GMT) 

Satellite 

09 

19:10:00 

NOAA17 

12 

18:01:00 

NOAA17 

12 

21:28:00 

NOAA18 

13 

19:18:00 

NOAA17 

14 

06:37:00 

NOAA17 

14 

18:55:00 

NOAA17 

14 

21:08:00 

NOAA18 

15 

18:32:00 

NOAA17 

15 

20:57:00 

NOAA18 

17 

19:26:00 

NOAA17 

17 

20:36:00 

NOAA18 

21 

10:11:00 

NOAA18 

21 

19:34:00 

NOAA17 

22 

19:11:00 

NOAA17 

22 

21:26:00 

NOAA18 

23 

18:48:00 

NOAA17 

23 

21:15:00 

NOAA18 

24 

18:25:00 

NOAA17 

24 

21:05:00 

NOAA18 

Table  III.  Satellite  Images  used  in  the  January,  2007  Case  Study 
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Date 

Time(GMT) 

Satellite 

17 

21:18:00 

NOAA18 

18 

19:21:00 

NOAA17 

18 

21:08:00 

NOAA18 

19 

18:58:00 

NOAA17 

19 

20:57:00 

NOAA18 

21 

18:12:00 

NOAA17 

21 

20:37:00 

NOAA18 

25 

21:36:00 

NOAA18 

26 

19:37:00 

NOAA17 

26 

21:26:00 

NOAA18 

27 

21:16:00 

NOAA18 

Table  IV.  Satellite  Images  used  in  the  October,  2006  Case  Study 
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